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Abstract 

A three year Cooperative Research Agreement 
between the Subsonic Aerodynamics Branch of 
the NASA Langley Research Center and the 
Virginia Polytechnic Institute and State Univer- 
sity (Va. Tech) has been completed. This doc- 
ument presents results from this three year en- 
deavour. 


Introduction 

The development of computational methods for 
predicting flight performance characteristics of 
vehicles has been the subject of intense re- 
search for many years. During this time, com- 
putational fluid dynamics methods for realis- 
tic aircraft configurations have been developed. 
However, these methods typically have focused 
on fixed wing aircraft operating in steady state 
flight conditions. Prediction methods for rotor- 
craft have lagged far behind fixed wing meth- 
ods due to the inherent unsteadiness and gen- 
eral complexity of rotorcraft flowfields. For ex- 
ample, in level, “steady” forward flight, rotor 
blades encounter unsteady flow due to the ro- 
tation of the blades. Even in hover, the flow- 
field is unsteady due to the asymmetry of the 
fuselage about the axis of rotation of the rotor 
blades. 

Many methods have been developed in the 
past for predicting aerodynamic interactions 
between a rotor and a fuselage with varying 
degrees of success. These methods are out- 
lined and referenced in the three Appendices to 
this document. In addition, the Appendices de- 
scribe in detail the method developed here. 

The present computational technique, which 
is described in the Appendices, was developed 


over a three year period. The following sections 
outline the research effort and the accomplish- 
ments made in each of the three years. 


Year 1 

The first year of the effort covers the pe- 
riod from November 1996 to November 1997. 
There were three primary goals set for this first 
period. These goals were geared toward feasi- 
bility studies, review of appropriate literature, 
and review of computer codes to determine the 
most promising path to take in the development 
of an efficient unsteady model. These goals 
were as follows: 

1. Review previous work leading to the 
pressure disk model used in the INS3D- 
UP incompressible Navier-Stokes code. 

2. Review the Generalized Dynamic Wake 
Theory (GDWT) to determine suitabil- 
ity for use as an unsteady pressure disk 
model coupled to a Navier-Stokes code. 

3. Examine possible candidate Navier- 
Stokes codes to determine suitability 
for time-averaged as well as time- 
accurate computations using a pressure 
disk model. 

All of these goals were accomplished in the 
first year. From the first item above, an un- 
derstanding of the current pressure disk model 
was obtained. For the second item above, a 
new computer code was written. This code, 
described in detail in the Appendices, was 
a developed specifically to compute the un- 
steady inflow and unsteady pressure jump to 
be used in conjunction with a Navier-Stokes 


computer code. Work on the third item con- 
sisted of examination of three Navier-Stokes 
codes. The first was INS3D, an incompress- 
ible code. While this is a versatile code, the 
convergence characteristics and CPU times of 
the code run in unsteady mode, while coupled 
to the GDWT, appeared to be unacceptable for 
the cases examined. The second code examined 
was CFL3D. Again, this is a versatile code and 
is widely used and well supported. Though the 
GDWT was never coupled to this code, discus- 
sions with the authors of CFL3D led to the con- 
clusions that the current overset grid capabili- 
ties and the time-accurate convergence charac- 
teristics of the code would deem it unacceptable 
for the purposes of this research. The third code 
examined was OVERFLOW. This code is well 
supported and has advanced overset grid capa- 
bilities, excellent time-accurate characteristics, 
and a low Mach number steady state capability 
which are appropriate for this research. Based 
on these rationales, OVERFLOW was chosen 
as the Navier-Stokes code for this research. 


Year 2 

The second year of the effort covers the pe- 
riod from November 1997 to November 1998. 
Again, there were three primary goals set for 
this second year period. While the first year 
goals were geared toward feasibility studies 
and literature review, the second year goals fo- 
cused on implementation of items from the first 
year. These goals were as follows: 

1. Begin coupling the GDWT to the OVER- 
FLOW code as a time-averaged pressure 
disk model. 

2. Implement the GDWT in the OVER- 
FLOW code as a time-accurate model. 


3. Investigate possible extensions to the 
GDWT/OVERFLOW model. 

All of the second year goals were completed 
in this year. Both the first and second items 
were completed and are detailed in the Ap- 
pendices. A number of extensions were ex- 
plored in this year as well. The paper in Ap- 
pendix I was written during this second year 
period. It describes the time-averaged and 
time-accurate model used for an isolated ro- 
tor. Time-averaged and time-accurate induced 
inflow predictions are shown to compare well 
with experimental data. Results from this pa- 
per were presented at the 54th Annual Forum 
of the American Helicopter Society in Wash- 
ington, D.C. in May 1998. 


Year 3 

The primary focus of the third year of the ef- 
fort was refining the model and making com- 
parisons between experimental data and pre- 
dictions. Appendix II is a Ph.D. dissertation 
that was written in this third year. This dis- 
sertation outlines in detail the method that was 
developed in the first two years. Predictions 
of time-averaged induced inflow, time-accurate 
induced inflow, steady surface pressures, time- 
accurate surface pressures for various configu- 
rations (an isolated fuselage, an isolated rotor, 
and a rotor/fuselage combination) are shown to 
compare well with experimental data. 


Summary 

In conclusion, a three year Cooperative Re- 
search Agreement has been completed. The 



goal of creating an efficient method to com- 
pute unsteady interactional effects between a 
helicopter rotor and fuselage has been accom- 
plished. The method is capable of predicting 
steady and unsteady induced inflow and steady 
and unsteady surface pressures on rotor and 
fuselage configurations. Currently, due to lim- 
itations that are introduced by almost all com- 
putational Navier-Stokes models, the predictive 
capability of the method is limited to moder- 
ate flight speeds. Note that this limitation is 
placed on the method by a particular compo- 
nent of the model and not by the overall method 
itself. Very low speed flight, including hover, 
will require further research into reducing or 
eliminating computational restrictions imposed 
by current Navier-Stokes codes. 
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Abstract 

An unsteady helicopter rotor model suitable for the 
study of rotor-fuselage interactional aerodynamics has 
been developed. The model couples a Generalized 
Dynamic Wake Theory (GDWT) code with a thin-layer 
Navier-Stokes code through an unsteady pressure jump 
boundary condition. The unsteady boundary condition 
models each rotor blade as a radially varying, time 
dependent pressure jump between adjacent grid planes. 
A non-rotating cylindrical grid is used for the isolated 
rotor case, and the unsteady boundary condition models 
the blades as a pressure jump traveling around the rotor 
disk. To demonstrate the coupling of the two codes, 
calculations are compared to experimental unsteady and 
time-averaged inflow data. 

Notation 


A fact 


e 0 

E add 

m, r 

[M] 

M 

[L c l [L s ] 


ratio of blade area in GDWT to area 
used in OVERFLOW, dimensionless 

local blade chord, dimensionless on R 
total energy per unit mass per unit vol- 
ume 

additional energy per unit mass per unit 
volume due to P , dimensionless 

harmonic number 
apparent mass matrix 
Mach number 

cosine and sine components of L- 
matrix 
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local load per unit span, dimensionless 
on pQ R 

polynomial number 

local pressure, dimensionless on P 

normalized Legendre functions 

radial location, dimensionless on R 
rotor radius, meters 

i-th component of perturbation velocity, 
dimensionless on SIR 
dimensionless time, dimensionless on 
Sl~ l 

freestream speed, dimensionless on 
SIR 

normal component of induced velocity, 
positive downward, dimensionless on 
SIR 

local mean normal component of 
induced velocity, positive downward, 
dimensionless on 

coordinate normal to rotor disk plane, 

dimensionless on R 

coefficients of the normal downwash 

series expansion 

effective angle of attack, rad 

mean coning angle, rad 
ratio of specific heats 
local blade pitch, rad 

mean and first harmonic components of 
local blade pitch, rad 
component of freestream normal to 
rotor disk, dimensionless on SIR posi- 
tive down 



^ advance ratio (component of freestream 

velocity parallel to rotor disk), dimen- 
sionless on QR 

v, T], \j/ ellipsoidal coordinates, dimensionless 

£ dimensionless coordinate along 

freestream line, positive pointing 
upstream from rotor 

P density, kg/m 3 

<X> pressure function, dimensionless on 

pQ 2 R 2 

qA qV acceleration and momentum compo- 

nents of pressure function, dimension- 
2 2 

less on p£2 R 

\j/ azimuth angle, rad/sec 

£2 rotor rotational speed, rad/sec 

( )* derivative with respect to t 

Introduction 

Unsteady interactional aerodynamics between a 
helicopter fuselage and rotor are not well understood. 
There exist a number of computational methodologies 
available to study these interactional effects. These 
range from superimposing isolated rotor and isolated 
fuselage effects to Navier-Stokes simulation of the 
entire problem [1-5]. Each of these methods have their 
respective advantages. For example, linearly 
superimposing the effects of an isolated rotor and 
isolated fuselage has the advantage of being a quick, 
simple method and works well when the nonlinear 
effects, such as flow separation, are negligible. On the 
other hand, Navier-Stokes simulations have the 
advantage of including all aspects of the flow field in a 
single calculation. Each of these methods also have 
disadvantages. For example, the superposition technique 
does not work well when the interactional effects are not 
truly linear, ( e.g . with flow separation), while full 
simulations may take many hours or even days of CPU 
time to execute. The hybrid method that is presented 
here draws from advantages of methods on either end of 
the spectrum. 

For typical shapes currently used for helicopter 
fuselages, flow separation is a common occurrence. This 
is typically the case even if one were to examine the 
isolated fuselage independent of the rotor system [6]. 
For isolated fuselage computations, a steady state 
solution is typically sought. Therefore, many of the 
solution acceleration techniques such as multigrid, grid 
sequencing, and preconditioning may be employed to 
obtain an economical solution to the isolated fuselage 


case. However, the inclusion of the rotor system to the 
computation adds even another level of complexity. If 
only the approximation to time-averaged rotor inflow 
effects are desired, several methods have been 
developed [3, 7, 8], If, however, viscous effects on a 
fuselage under the influence of unsteady rotor inflow 
effects are needed when separation is expected, only 
Navier-Stokes simulations of the entire problem 
including time-accurate representation of the moving 
blade geometry are presently available. For example, 
Meakin [5] used the unsteady, thin-layer Navier-Stokes 
equations to compute the flow field around a complete 
tiltrotor aircraft, including three-bladed rotors and 
moving overset grid systems for a hypothetical flight 
condition. Ahmad, et al, [9] also used the thin-layer 
Navier-Stokes equations to compute the unsteady blade 
pressures and flowfield on an isolated two-bladed rotor 
using moving overset embedded grids. In [6, 9] a first 
order time accurate scheme was used which necessitated 
the use of very small time steps. Even though solutions 
in [9] show a converged (periodic) solution in three to 
five rotor revolutions, 1 152 time steps were required per 
revolution at a cost of 45 Cray C-90 hours. 

The current work is designed to examine techniques 
to reduce the computational times required to provide 
unsteady calculations that include a rotor model while 
maintaining the capability to calculate viscous flows on 
a fuselage. It is not to replace the full simulations or the 
simpler analyses, but rather to provided an additional 
tool to the researcher or designer. 

Method 

The current method employs a hybrid approach 
which allows a quick assessment of the time-averaged 
and time-accurate rotor-fuselage interactional 
aerodynamics. A simple, fast-running Generalized 
Dynamic Wake Theory [10, 11] method is used to trim 
an isolated rotor and generate a time dependent pressure 
distribution to represent the rotor blades in the plane of 
the rotor disk. This set of unsteady, periodic pressures is 
used in a version of the thin-layer Navier-Stokes code, 
OVERFLOW [12], which was modified under the 
current effort to include an unsteady pressure jump 
boundary condition between adjacent planes in a non- 
rotating cylindrical grid for the isolated rotor. 

GDWT 

Details of the GDWT can be found in [10, 11]; 
however, for completeness, an outline of the theory is 
given here. In addition, the current method of time 
integration of the GDWT equations is presented as is the 
current method of trimming the rotor to specified thrust, 
roll, and pitch moments. 
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The GDWT [10, 11] begins with the mass 
conservation equation and the incompressible, 
linearized Euler equations: 

?U = 0 (1) 

W 

where repeated subscripts imply that the summation 
convention is being used and commas imply that 

derivatives are being taken. The pressure function, <£, 

A V 

can be split into two components, <1> and 

corresponding respectively to the two terms on the left 
hand side of equation (2), each of which satisfies 
Laplace’s equation. The boundary conditions are that 
the pressure function matches the known blade loading 
at the rotor and that it equals zero at infinity. 

When written in ellipsoidal coordinates, the 
solution to Laplace’s equation for the pressure potential 
on a circular planform, along with the conditions that 
the pressure perturbation equals zero at infinity and that 
the pressure potential is zero at the edge of the planform, 
can be written in closed form as follows: 


® = £ P>)T n {i n) 

m = 0 /i = m + 1, m + 3, ... 

[t^ c cos(m\j/) + T™sin(m\j/)] (3) 
where v , r \ , and \j ) are dimensionless ellipsoidal 
coordinates, P and Q are normalized Legendre 


functions, t™ c and are the unknown coefficients of 
the cosine and sine functions and are themselves 
functions of the dimensionless time variable, t . Since 

<X> A and each satisfy Laplace’s equation, equation 


(3) can be applied to and <t> v separately. Examining 
only the normal component of the perturbation velocity 
and integrating along a streamline from a point on the 
rotor disk back to infinity, leads to the following 
equations for the normal component of perturbation 
velocity: 


1 r 3<s> * 

w = --Jo TT* 


(4) 


dz 


(5) 


It] = 0 


where ( )* denotes a derivative with respect to 
dimensionless time and z is the direction normal to the 
rotor disk. Expanding the w component of velocity in 
an expansion similar to that used for O and 

superimposing the <$ A and <t> V components results in 
the following set of first order differential equations for 
the unknown, time-dependent coefficients (“states” of 


the model), a and (3, of the w velocity component 
expansion: 
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Closed form expressions for [A/] , [L c ] , [L 5 ] , x" c , 


and x™ 5 are given in [10, 1 1]. As described in references 

10 and 1 1, the t functions depend on the loading model 
chosen by the user. The following loading model is 
used: 


[w,-X + P 0 ncosy-y 
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0 = 0 O + 0 c cosy + e^iny 
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(9) 
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V 1 -M 2 

This loading model includes angle of attack effects due 
to (1) the freestream component of inflow normal to the 
disk, (2) the induced inflow feedback from the GDWT, 
(3) the mean coning angle of the blades, (4) the pitch 
rate, (5) the tangential velocity variations due rotor 
rotation, and (6) the increase in lift with Mach number. 
In equation (10), a lift curve slope of 2k has been 
assumed. A very simple stall model is also included in 
the method whereby the effective angle of attack is not 
allowed to increase above 10 degrees. However, it has 
been found that, for the cases presented here, the results 
from this model are not sensitive to the 10 degree 
maximum angle of attack. 

Also described in these references is a technique to 
generalize equations (6) and (7) to include effects of the 
mean induced inflow on the wake skew angle. It is this 
non-linear version of the GDWT that is used in the 
results presented in this paper. 

Equations (6) and (7) form a set of first order 
ordinary differential equations that can be solved by 
various techniques. Here, a Jameson-style Runge-Kutta 
integration [13] is used to solve for a and (3 , and thus 
for w , at each successive time step. Time integration is 
carried out over a specified number of rotor revolutions. 
That leads to the complete loading and downwash 
solution on the rotor disk. This loading information is in 
turn used in an outer trim loop for the isolated rotor. 

The rotor trim condition is determined using a 
modified Newton-Raphson technique [14]. In this 
technique, the rotor state is first determined based on a 
specified collective, lateral, and longitudinal pitch 
setting. These controls are then each perturbed 



independently to form a “derivative matrix” which is 
used to determine the successive collective, lateral, and 
longitudinal pitch settings for the rotor based on the 
changes noted in the rotor thrust, roll, and pitch 
moments due to the control perturbations. These 
successive determinations of the pitch settings are 
continued until the rotor is trimmed to within a specified 
root-mean-square change in the thrust, roll, and pitch 
moments. This trim technique is classified as a modified 
Newton-Raphson technique since the derivative matrix 
is held fixed throughout the trim procedure rather than 
being recalculated at each step of the iteration as it 
would be done in a classic Newton-Raphson method. 

Once the rotor has been trimmed, the unsteady 
loading and unsteady downwash are known for the 
entire rotor disk. The unsteady downwash has a 
frequency content consistent with the specified number 
of harmonics and shape functions chosen by the user in 
the GDWT method. Four harmonics provide good 
correlation between predicted and measured downwash 
when using a loading model that assumes the blades are 
modeled as a “pressure delta function” or pressure spike 
traveling around the rotor azimuth. This is the method 
used in this paper. In addition, since the non-linear 
version of the GDWT is used, the calculated loading is 
coupled to the calculated downwash. Thus, the loading, 
which initially was two-dimensional, has now been 
corrected to account for effects up to a frequency of four 
per revolution. It is this discrete, unsteady rotor disk 
pressure distribution which is to be used as a boundary 
condition in the thin-layer Navier-Stokes code. 

OVERFLOW 

In this section, the thin-layer Navier Stokes code, 
OVERFLOW, and the method used to couple 
OVERFLOW to the GDWT code is discussed. 

For purposes of this study, the thin-layer Navier- 
Stokes code, OVERFLOW (version 1.7v), was chosen 
as a baseline code based on its robustness, its 
convergence times for unsteady cases, and its low Mach 
number capabilities. For this study, this version of 
OVERFLOW was modified to include an unsteady 
pressure jump boundary condition. This boundary 
condition is applied between two user specified planes 
which are separated by a user specified “iblanked” 
plane. (As is common practice, an “iblanked” region is 
one which is assumed to be outside the computational 
domain and in which no computations are performed.) 
The new boundary condition is implemented as an 
additional energy term which is added to the quantity 
pe 0 in the following form: 


where P g is the dimensionless pressure determined by 

the local load per unit span (from the GDWT code) 
divided by the local blade segment chord at the current 
grid point, y is the ratio of specific heats, and Aj- acl is 

the ratio of the total blade area in the GDWT to the total 
area on which the pressure is applied in OVERFLOW in 
order to maintain the same overall thrust between the 
two methods. For the two specified planes where the 
boundary condition is applied, one half of the additional 
energy term is placed on each plane. The code is then 
executed in the time-accurate mode until a periodic 
inflow solution is obtained. 

Code Coupling 

In the current effort, the GDWT code and 
OVERFLOW are used together in a “loose” coupling. 
That is, the GDWT code and OVERFLOW are both 
stand-alone codes that are executed sequentially. Since 
the GDWT code calculates the isolated rotor trim 
loading, a trim calculation is not needed in 
OVERFLOW. Thus, OVERFLOW is required only to 
execute enough time steps to obtain a periodic solution. 

The common link between the two codes is that the 
same rotor system pressure distribution is used in each; 
the GDWT calculates a pressure distribution that is 
subsequently used in OVERFLOW. However, there is a 
difference in how the pressure distributions are used in 
each code. The GDWT, as applied here, assumes a that 
the loading is a delta function at the mid-chord of the 
blade, whereas OVERFLOW assumes that the loading is 
applied on grid line which has associated with it an area 
“wedge” in the cylindrical grid used for the rotor. In 
addition, the GDWT uses the linearized Euler equations 
and OVERFLOW uses the Navier-Stokes equations. 
Therefore, the common pressure distribution when used 
in each code will produce slightly different inflow 
characteristics due to the slightly different physical 
assumptions. Since the GDWT does not include a 
fuselage model, it will be necessary to iterate between 
the two codes when a fuselage is introduced into the 
OVERFLOW solution. An iterative capability has been 
implemented in the codes, but is not necessary in the 
present work since only the isolated rotor case is shown. 

In order to include fuselage effects in the GDWT, 
an iterative technique has been developed in which the 
inflow from the GDWT and OVERFLOW are compared 
on a similar frequency basis. Since the GDWT provides 
inflow information only up to user specified frequency, 
the inflow calculated by OVERFLOW is filtered to 
match the GDWT frequency for purposes of the 
comparison. A difference is taken between the GDWT 
inflow and the filtered OVERFLOW inflow. This 
difference is then used in the GDWT code as an “inflow 


correction” to account for differences in the two codes 
and assumptions made therein. The GDWT code is re- 
executed, and the entire procedure may be repeated as 
needed. 

Experimental Data 

All of the experimental data used in this paper are 
from laser velocimeter inflow measurements obtained at 
the NASA Langley 14- by 22-Foot Tunnel [15]. Two 
different blade planforms were used in the test, one 
tapered and one rectangular. Only the data for a rotor 
with a 4-bladed, tapered planform is presented in this 
paper. This rotor has constant chord (3.15 inches) over 
the inboard 75 percent of the span and an approximately 
3-to-l taper ratio over the outer 25 percent of the span 
such that the tip chord is approximately one inch. The 
blades have a linear twist of -13 degrees, a root cutout at 
25 percent span, a solidity of 0.0977 and radius of 32.5 
inches. The thrust coefficient was 0.0065. Inflow data 
used here were taken in a plane approximately 2.6 
inches above the tip path plane of the rotor for all 
azimuth stations 30 degrees apart and at a number of 
radial stations. The flight condition chosen for this 
comparison has a 3 degree nose down shaft tilt and an 
advance ratio of 0.23. Data samples were processed at a 
resolution of 128 samples per revolution or about every 
2.8 degrees of blade travel. 

Results 

GDWT Time-Averaged Results 

In this section, results will be presented from the 
GDWT code before the coupling with OVERFLOW. 

Figure 1 shows a contour plot comparison between 
the time-averaged, measured inflow ratio and the 
GDWT code prediction for the baseline case. For the 
GDWT calculations, 4 harmonics and 15 states were 
used. It can be seen that many of the major time- 
averaged flow features are captured, as discussed in 
[16]. An upwash region can be see along the front of the 
disk in both the measurement and in the prediction. The 
prediction matches the location of the zero induced 
inflow ratio line well, but the predicted gradients of the 
inflow, especially on the advancing and retreating sides, 
are higher than the gradients of the inflow in the 
measured data. Other features such as the strong 
downwash region in the first quadrant are predicted 
well. Some of the discrepancies in the inflow 
distribution can be attributed to the fact that there was a 
fuselage body and hub present when the experiment was 
conducted [16], but the GDWT in this case includes no 
fuselage or hub effects. 


GDWT Time-Accurate Results 

Figure 2 shows a comparison of the measured and 
predicted unsteady inflow ratio with mean inflow ratios 
removed for nine experimental measurement locations. 
The measurements are located above the rotor plane 
along 3 radial lines at azimuth angles of 30, 180, and 
300 degrees and at 3 spanwise locations of r/R=0.60, 
0.78, 0.90. These measurement locations are marked 
with the letters A-J (excluding the letter I) on figure 1. In 
all of the plots in figure 2, the measured and predicted 
data show a general four per revolution oscillation in the 
inflow as is expected for a four bladed rotor. The 
locations D-F (r/R=0.78) in this figure show that 
predictions of unsteady inflow peak-to-peak magnitudes 
and phases are predicted well. For all of the plots in this 
figure, good phase predictions are displayed, but peak- 
to-peak predictions are not as good. The measured data 
generally show sharp blade passage pulses at the 
beginning of each oscillation cycle. The predictions 
using the GDWT in this manner are not expected to 
predict these sharp blade passage pulses, since only the 
first four harmonics of inflow are being computed (i.e. f 
only 15 states in the GDWT model are used). However, 
the discrete loading distribution used in the GDWT to 
calculate the inflow is not as strongly affected by the 
number of harmonics used. It is these loading values 
that are to be used in OVERFLOW. It has been shown 
that using many more states can improve the 
correlations of the time- accurate inflow [17]; but the 
corresponding improvement in the loading distribution 
is much smaller and hence much less necessary. The 
results presented in figures 1-2 give a general confidence 
that the GDWT code is working correctly. 

OVERFLOW Time-Averaged Results 

With the loading predictions from the GDWT code, 
OVERFLOW is executed in two stages. First, to 
minimize the number of time steps required to converge 
to a periodic solution, a mean flow is established by 
using a steady-state computation with a constant 
pressure jump equivalent to an equal distribution of 
thrust on the rotor. Then, the case is restarted in time- 
accurate mode using the GDWT determined unsteady 
pressure jump boundary condition. For the cases 
presented here, the grid topology consists of a single 
stretched cylindrical grid which extends 1.5 rotor radii 
above and below, 2.5 rotor radii ahead and to either side 
of, and 4.5 rotor radii behind the rotor (see figure 3). A 
grid size of 115x129x56 was used in the vertical, 
azimuthal, and radial directions, respectively. The rotor 
portion of the grid is in the vertical center of the grid and 
occupies a grid of 129x30 (minus the inner 6 radial 
stations, to account for the blade root cutout). A 
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freestream/characteristic (Reimann invariant) boundary 
condition was used at the outer boundaries; a periodic 
boundary condition was used at the 0 degrees 

location; an axis boundary condition was used along the 
rotor axis; the new unsteady boundary condition was 
used along two planes in the vertical center of the grid 
(separated by a single, iblanked plane of the same size 
as the rotor planes) and extends from the root cutout 
location to the rotor edge. For all of the cases examined 
here, 1408 time steps (a multiple of the number of time 
steps per revolution) were sufficient to obtain a 
converged steady state solution. Here, convergence was 
determined by a two order of magnitude drop in the L2- 
norm of the residual. Though typically a larger residual 
drop considered crucial to a converged solution, a two 
orders of magnitude drop in residual is sufficient here 
due to the small values of the initial residual. This 
steady state calculation was then followed by two 
revolutions of unsteady calculations, which was 
sufficient to obtain periodicity in the solution. 

The computations presented use the following 
methods in OVERFLOW: (1) central difference 
calculations of the right-hand side of the equations, (2) 
3-factor diagonal left-hand side inversion, (3) a matrix 
dissipation scheme, and (4) Newton sub-iterations. 

Figure 4 shows a contour plot comparison between 
the time- averaged, measured inflow ratio and the 
OVERFLOW prediction which has been filtered to 
contain the same frequency content as the GDWT 
prediction presented in figure 1. For purposes of 
comparisons later in the paper, this case is referred to as 
the “baseline case”. Since the highest frequency 
available in the measured time-averaged contour plot is 
six per revolution (due to Nyquist cut-off 
considerations), this filtering also places the 

OVERFLOW prediction on a frequency basis similar to 
those measured data. Generally, the major time- 
averaged features seen in figure 1 are also seen in the 
figure 4. It can be seen in figure 4 that the location of the 
zero induced inflow line is not as well predicted, but the 
overall levels are well predicted. Also, in contrast to the 
results presented in figure 1, the gradients of inflow 
(spacings between the contour lines) are predicted well. 
As with the GDWT time-averaged results presented 
earlier, some of the discrepancies in the flow field can be 
attributed to the fact that the experimental data includes 
the effects of a fuselage and hub. An example of these 
discrepancies can be seen on the retreating side of the 
rotor disk near the center of the contour plot. At that 
location in the predicted results, a region of near zero 
inflow can be seen which is not present in the measured 
data. Since the predictions do not use a hub or hub wake 
model, there are no hub blockage effects in the flowfield 
and the flow is free to travel through the center of the 


rotor disk where there is no blade loading. Since it is the 
rotor loading and hub/fuselage blockage effects that 
actually cause the induced inflow flowfield, it can be 
deduced that the near zero inflow region near the center 
of the rotor disk is at least partially caused by the low 
loading and lack of hub model in the hub region. 

OVERFLOW Time-Accurate Results 

Figure 5 shows a comparison of the measured 
unsteady inflow and the predicted unsteady inflow. In 
these plots, the sharp blade-passage pulses are seen at 
the beginning of each oscillation cycle, and the 
magnitudes and phases of these pulses are predicted 
well. 

Sensitivities 

To show the sensitivity of the above results* to 
particular parameters related to the solution procedure, 
several of the relevant parameters pertaining to the 
solution procedure are explored. They include the 
number of Newton sub-iterations, the use of viscous 
terms, the time step (azimuthal resolution) used, and the 
outer boundary condition used. In addition, the amount 
of CPU time required for each case is presented. 

Newton Sub-iterations 

In the case presented above, six Newton sub- 
iterations were used at each time step in the unsteady 
portion of the calculations. To demonstrate that six sub- 
iterations is sufficient, the baseline case was re-executed 
with ten sub-iterations (instead of six) starting from the 
same steady state conditions. Figure 6 shows the 
normalized L2-norm of the conservative variables 
versus Newton sub-iteration number for the last time 
step in the case. It can be seen that there is little 
advantage to using more than six Newton sub-iterations. 
This is further demonstrated by figures 7 and 8, which 
show that there is virtually no difference in the time- 
averaged or time-accurate inflows. Thus, the use of six 
Newton sub-iterations is adequate for the current case. 

Viscous Terms 

Since the case presented here is an isolated rotor 
with no surface in the flowfield, viscous terms in the 
flow solver were not used. To show that these terms are 
not dominant in this particular case, the above baseline 
case was re-executed with all of the thin-layer viscous 
terms activated in the flow solver. Figures 9-10 
demonstrate that these terms are not a dominant 
influence in the isolated rotor case. However, these 
viscous terms will be required when a fuselage or other 
solid surface is introduced into the solution. 



Table 1: 


Time Step 

All of the above cases used 128 time steps per 
revolution. To test the adequacy of this time step (and 
azimuthal resolution), the entire case were re-executed 
with a doubled time step, i.e., with 64 azimuth steps per 
revolution. The resulting time averaged and time 
accurate calculations are shown in figures 11-12. Both 
of these figures show that 64 time steps per revolution is 
adequate to capture the gross features of the time 
averaged and time accurate inflow. Due to the 
similarities between the 64 time step per revolution case 
and the original 128 time step per revolution case, it 
appears that the 128 time step per revolution case is 
quite adequate for these predictions. 

Outer Boundary 

In all of the above cases, the outer boundary is 
relatively close in proximity to the rotor disk. This raises 
a question as to the effect of the outer boundary on the 
solution. To demonstrate the effect of the outer 
boundary on the solution, the boundary condition there 
was changed from a Reimann invariant (freestream/ 
characteristic) condition to simply a freestream 
condition and the baseline case was re-executed. Figures 
13-14 show the results of the computations with the 
freestream outer boundary condition. Figure 13 shows a 
good match with the baseline time averaged plot. In 
figure 14, it can be seen that there is a slight blade-to- 
blade difference at the beginning of the revolution which 
disappears by the end of the revolution. This is evidence 
that, formally, one more revolution should be calculated, 
however, for the remaining blade passages, it can be 
seen that the gross effects of the outer boundary 
condition on the unsteady inflow are small. This implies 
that the solution is relatively insensitive to what is 
happening at the outer boundary. 

CPU times 

All of the steady and unsteady calculations using 


OVERFLOW were run on a Cray C-90. Table 1 shows a 
comparison of the CPU times for each of the cases 
presented above. As can be seen in the table, the CPU 
time required is relatively inexpensive. 

Table 1: 


Case 

CPU (hr: min) 

Baseline 

4:16 

Newton Sub-iteration 

5:39 

Viscous Terms 

5:01 


Case 

CPU (hr: min) 

Time Step 

1:38 

Outer Boundary 

4:06 


Concluding Remarks 

(1) An unsteady helicopter rotor model using an 
unsteady pressure jump boundary condition has been 
developed and included in the thin-layer Navier-Stokes, 
overset grid code, OVERFLOW. 

(2) Even though isolated rotor calculations are 
presented here for the purposes of demonstrating the 
model, excellent agreement is shown between time 
averaged and unsteady measured and predicted inflow 
ratios at a plane above the rotor disk, and CPU 
requirements for the isolated rotor method are found to 
be reasonable. 

(3) It has been demonstrated that suitable 
parameters were chosen with regard to the number of 
Newton sub-iterations used, the viscous terms used, the 
time step used, and the outer boundary conditions used. 
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Figure 2: Measured and GDWT Predicted Time Histories of Induced Inflow for Azimuth Angles of 30, 180, and 
300 degrees at Three Radial Stations of r/R = 0.60, 0.78, 0.90 


Figure 3: Schematic of Isolated Rotor Grid Topology as Viewed from Behind and Above Rotor on the Retreating 
Side 
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Figure 5: Measured and OVERFLOW-GDWT Predicted (Baseline) Time Histories of Induced Inflow for Azimuth 
Angles of 30, 180, and 300 degrees at Three Radial Stations of r/R = 0.60, 0.78, 0.90 
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Figure 6: Normalized L-2 norm of Conservative Variables vs. Newton Sub-iteration Number 







Inflow Ratio - Mean Inflow Ratio - Mean Inflow Ratio - Mean 





Ref. Blade Location 


r/R = 0.60 
Location B 



r/R = 0.78 



r/R = 0.90 
Location H 



Ref. Blade Location 


Baseline 
Baseline with 10 
Newton Subiterations 




Location J 



Ref. Blade Location 


Figure 8: Time Histories of Induced Inflow for Baseline vs. Baseline with 10 Newton Sub-iteration for Azimuth 
Angles of 30, 180, and 300 degrees at Three Radial Stations of r/R = 0.60, 0.78, 0.90 
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Figure 10: Time Histories of Induced Inflow for Baseline vs. Baseline with Viscous Terms Included for Azimuth 
Angles of 30, 180, and 300 degrees at Three Radial Stations of r/R = 0.60, 0.78, 0.90 
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Figure 12: Time Histories of Induced Inflow for Baseline vs. Baseline with Doubled Time Step (Halved Azi- 
muthal Resolution) for Azimuth Angles of 30, 180, and 300 degrees at Three Radial Stations of 
r/R = 0.60, 0.78, 0.90 
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ROTOR/FUSELAGE UNSTEADY 
INTERACTIONAL AERODYNAMICS: 
A NEW COMPUTATIONAL MODEL 


David Douglas Boyd, Jr. 


ABSTRACT 


A new unsteady rotor/fuselage interactional aerodynamics model has been developed. This 
model loosely couples a Generalized Dynamic Wake Theory (GDWT) to a Navier-Stokes solution 
procedure. This coupling is achieved using a newly developed unsteady pressure jump bound- 
ary condition in the Navier-Stokes model. The new unsteady pressure jump boundary condition 
models each rotor blade as a moving pressure jump which travels around the rotor azimuth and 
is applied between two adjacent planes in a cylindrical, non-rotating grid. Comparisons are made 
between predictions using this new model and experiments for an isolated rotor and for a coupled 
rotor/fuselage configuration. 
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English symbols 

000 freestream speed of sound [m/sec] 

a ™ series expansion coefficients [see Appendix A] 

A(r) local area ratio between blade and computational cell 

1 the imaginary number, y/—l 

c local blade chord normalized by R 

c coefficient function [see Appendix A] 

Cf,C m force and moment coefficients 
C p pressure coefficient 

d p modified pressure coefficient 

dp unsteady component of modified pressure coefficient 

Cj thrust coefficient 

Ct mean thrust coefficient 

Q, roll moment coefficient 

Cl mean roll moment coefficient 

Cm pitch moment coefficient 

Cm mean pitch moment coefficient 

eo stagnation energy per unit mass [m 2 /sec 2 ] 

F, G, H inviscid fluxes (see equations (4.3) to (4.5)) 

F V ,G V ,H V viscous fluxes (see equations (4.3) to (4.5)) 
h hinge offset location normalized by R 
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coefficient function [see Appendix A] 
quasi-steady inflow matrix, cosine component 
quasi-steady inflow matrix, sine component 
sectional blade lift [N/m] 
local Mach number 
freestream Mach number 

number of azimuthal time steps per rotor revolution 
pressure [N/m 2 ] 
freestream pressure [N/m 2 ] 

normalized, associated Legendre functions of the 1st kind 

normalized, associated Legendre functions of the 2nd kind 

i th component of perturbation velocity [m/sec] 

cartesian components of heat flux [J/(m 2 sec)] 

vector of conservative variables = [p, p u, pv, pw, peo] r 

i ,h conservative variable 

radial location on disk, normalized by R 

radial location on disk [m] 

rotor radius [m] 

gas constant [J/(kg K)] 

time [sec] 

temperature [AH 

nondimensional time, normalized by £2 
freestream velocity, normalized by Q.R 
induced inflow velocity [m/sec] 
elements of mass flow matrix 
cartesian velocities [m/sec] 

local normal component of velocity at rotor disk, normalized by SIR 
local normal component of w 
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Greek symbols 


a eff 

a* 

Po 

Plc»Plj 

Po 

Vj&J 

Y 
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AP 

8 

X 
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^i^mean 

A(r,\|r) 

<X> 

<(r) 

V,T1,V 

V 
Vi 


Pi 

Hi } mean 

0 

0rw(^) 

00 

0(T) 9s 

01 


effective angle of attack [rad] 

rotor shaft tilt angle [rad] 

blade mean coning angle [rad] 

first harmonics of rotor flapping [rad] 

blade mean coning angle [rad] 

induced inflow coefficients 

ratio of specific heats 

Dirac delta function [see Appendix A] 

pressure jump [N/m 2 ] 

user specified tolerance 

freestream component normal to disk, normalized by Q.R, positive down 
induced inflow ratio = V;/ ( Q.R ) 
time averaged induced inflow ratio 
mean component of A,- 

generic time averaged function [see Appendix A] 
pressure function, normalized by p Q?R 2 
radial expansion shape function 
dimensionless ellipsoidal coordinates 
ijlr on rotor disk 

phase of first harmonic of blade pitch [rad] 
local air density [kg/m 3 ] 
freestream air density [kg/m 3 ] 
rotor advance ratio = V^/QR 

induced inflow in fu direction, normalized by QR, positive downstream 

mean component of /j,- 

local blade pitch [rad] 

built-in blade twist function [rad] 

blade collective pitch [rad] 

cosine and sine components of blade pitch [rad] 

magnitude of first blade pitch harmonic [rad] 
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cosine part of pressure expansion of m th harmonic and n lh polynomial 
sine part of pressure expansion of m th harmonic and n th polynomial 
stress tensor [m/sec] 

rigid flap natural frequency [cycles/revolution] 
rotor rotational speed [rad/sec] 


Subscripts 


additional component 
i th component 

derivative with respect to i ,h direction 
polynomial number 

coordinate along freestream line, positive pointing upstream 
freestream quantities 


Superscripts 


harmonic number 
unsteady component 



Operators & Acronyms 
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A(-) operator used to denote an incremental value of (•) 

^F(-) filter operator 

(•) ! ! double factorial [see Appendix A] 

* 

(•) derivative with respect to t 

(•) reference quantities 

GDWT Generalized Dynamic Wake Theory 

RIM rotor loading model 

RFFM rotor/fuselage flowfield model 
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Chapter 1 


Introduction 


1.1 Motivation 

It is well known that rotary wing aircraft aerodynamics are complicated. Unlike fixed wing aircraft, 
on which a steady-state flight condition typically implies steady-state aerodynamics, a rotary wing 
aircraft experiences a significant unsteady aerodynamic environment in all flight conditions, even 
in level, unaccelerated flight, due to the presence of the rotating wings (rotors). This aerodynamic 
environment includes the aerodynamic interactions, which are inherently unsteady and complex, 
between the rotor(s) and the fuselage. One example of the complexity associated with these inter- 
actions is the problem of flow separation phenomena. Whereas fixed wing aircraft typically have 
little significant flow separation in steady-state flight due to their streamlined fuselage shapes, ro- 
tary wing aircraft typically have blunt aft regions that are conducive to flow separation. Even in 
hover, the flow induced by the rotor(s) impinging on the fuselage tends to separate on the underside 
of the fuselage which is often a blunt surface. Sheridan and Smith [1] discuss many other examples 
and categories of rotorcraft interactional aerodynamics. They also state in their conclusions that. 

“ . . it will be necessary to develop tractable theories and analytical methods to account 
for all these phenomena. Interactional aerodynamics of the airframe is not as neatly 
packaged as rotor aerodynamics. Many of the interactions involve viscous processes, 
and in some aspect semi-empirical techniques may always be needed. But a start must 
be made in developing the required mathematical models so that we can cope with 
these problems adequately in the vehicle design phase.” 
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In short, the above passage calls for the development of methods to address the coupled ro- 
tor/fuselage interactional aerodynamic effects, whereas previously, the rotor effects and fuselage 
effects were treated in isolation. 

In the years since Sheridan and Smith’s paper [I], many types of analyses have been developed 
and used in the prediction of the unsteady interactional aerodynamic characteristics of rotorcraft. 
Figure 1.1 is a graphic that depicts and categorizes several of the major methodologies. 

In the area of relatively low computational expense and complexity (see figure 1.1), singularity 
methods have been used. These methods typically use singularities, such as a lifting line to repre- 
sent the rotor blade, a system of vortices to represent the wake, and source, doublet, and/or vortex 
panels to model the fuselage. In a relatively inexpensive and computationally efficient manner, 
these methods are able to capture low order effects on each component due to the other compo- 
nent, such as the mean downwash on fuselage due to the rotor or the mean inflow at the rotor disk 
due to the fuselage. But, since the fuselage is typically modeled using a panel method, calcula- 
tion of some interactional aerodynamic effects, such as flow separation due to rotor downwash, is 
difficult. In cases where viscous effects are predominant, the viscous flow effects must be either 
ignored, specified a priori , or determined by coupling the method to a boundary layer type model. 

On the other end of the computational expense and complexity scale, at relatively high computa- 
tional expense and complexity, are the methods involving computational fluid dynamics (CFD), in 
particular, Navier-Stokes methods. These types of methods have been used to calculate the entire 
flow field of the complete rotorcraft configuration, all in one computation. Even though these com- 
putational methods are theoretically able to capture all of the interactional aerodynamic couplings 
between the rotor(s) and fuselage, their computational expense is prohibitive for routine use. 

There is a lack of methods available in the literature which fall between the singularity methods 
and the Navier-Stokes CFD methods for studying unsteady interactional aerodynamics of rotor- 
craft. The current research is motivated by the lack of available hybrid methods which are compu- 
tationally efficient, yet are able to capture primary interactional aerodynamic effects between the 
rotor and fuselage. Figure 1.1 shows, with a dotted ellipse labeled “Hybrid Methods , where the 
current work falls on the computational expense and complexity scale. 
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1.2 Literature Review 

As discussed previously, unsteady rotor/fuselage interactional aerodynamics generally fall into 
three categories: (1) singularity methods, (2) CFD methods, and (3) hybrid methods. In the fol- 
lowing sections, a brief review of each is given. 

1.2,1 Singularity Methods 

Singularity methods are typically characterized by the use of a source, doublet, and/or vortex panel 
representation of the fuselage, a lifting line or lifting surface representation of the rotor, and a vor- 
tex lattice model representation of the rotor wake. For rotorcraft analyses, these methods are used 
to compute the flowfield of the complete vehicle. Johnson [2] provides an extensive discussion of 
singularity methods used for rotorcraft analyses up through the year 1986. Since that time, other 
singularity methods have been developed as well. Egolf and Lorber [3] used a source panel de- 
scription of the fuselage, a lifting line blade model, and a prescribed vortex wake description of 
the rotor/fuselage system to model the unsteady rotor/fuselage flowfield. The prescribed vortex 
wake was prevented from cutting through the fuselage by displacing, in an a priori manner, the 
segments of the vortices that would otherwise have been inside the fuselage. No attempts were 
made to model the wake of the fuselage or the flow separation from the fuselage. Only limited 
comparisons were made to experimental data. Mavris, et al. [4] used a doublet representation of 
the fuselage, a lifting line blade model, and a free vortex wake description of the rotor/fuselage 
system. No modeling was used for the fuselage wake or fuselage flow separation. Also, vortex 
wake filaments that are inside the fuselage were excluded from the computations. Comparisons 
with experimental pressures show good agreement along the top of the fuselage, but agreement 
degrades on the sides of the fuselage. Mavris, et al. [4] attribute these discrepancies to flow sepa- 
ration on the fuselage and to inadequate vortex-surface interaction predictions. Berry [5] combined 
a fuselage source panel representation, a source-dipole representation for the rotor, and a distorting 
vortex-lattice representation of the rotor wake to model the rotor/fuselage system. Comparisons are 
made to measured time averaged and unsteady inflow velocities; no fuselage pressure comparisons 
are made which include the rotor influence. Quackenbush, et al. [6] used a source/doublet de- 
scription of the fuselage, a vortex-lattice model for the rotor blades and a novel “Constant Vorticity 
Contour (CVC)” free wake model; fuselage flow separation was not modeled. Close surface/vortex 
interactions were modeled using selective remeshing of the curved vortex elements and using an 
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“Analytical/Numerical Matching (ANM)” scheme. Computational efficiency was improved by us- 
ing “fast vortex methods” for wake-on-wake and wake-on-body computations. Generally good 
agreement with measured results are demonstrated for time averaged induced velocities above the 
rotor disk and for time averaged and unsteady pressures on the top centerline of the fuselage. 
Crouse [7] used a source panel description of the fuselage, a lifting line blade model, and a free tip 
vortex wake model without an inboard wake model to represent the rotor/fuselage system. Vortex 
wake elements that cross the fuselage surface are handled by splitting them into smaller segments, 
and shifting the collocation points of these smaller segments such that they are at a specified min- 
imum distance from the fuselage surface. This method is similar in concept to that used by Egolf 
and Lorber [3] as discussed above. Good comparisons of unsteady pressures were shown on the 
top centerline of the tail boom of the fuselage. Boyd, et ol. [8] included the open-loop effects of 
a fuselage, represented by a non-lifting fuselage source panel method, on the rotorcraft trim in a 
comprehensive rotorcraft analysis. This effect was implemented in the comprehensive analysis as 
an additional rotor inflow distribution plus an additional rotor wake distortion due to the presence 
of the fuselage. Effects of the rotor on the fuselage were not modeled. Though some computations 
have proved successful and can be computationally efficient, all of these singularity methods suffer 
from the inability to predict some of the rotor/fuselage interactional effects. For example, methods 
that use a source panel description of the fuselage do not have the capability to determine the lift 
or the lift change on a fuselage due to the rotor. Also, quantities such as flow separation and drag 
must either be ignored, must be specified a priori , or must be determined by coupling the method 
to a boundary layer model. 


1.2.2 CFD Methods 

In recent years, unsteady calculations on complete rotorcraft configurations using CFD methods 
have become possible. In addition, there are several degrees of complexity that can be modeled 
with CFD. For rotorcraft applications, methods have been developed to solve the full potential 
equation, the Euler equations, and the Navier-Stokes equations. 

Chen and Bridgeman [9] coupled the three dimensional boundary layer equations to the full 
potential equations in a blade-fixed coordinate system for an isolated rotor. The three dimensional 
boundary layer equations assumed that the surface curvature effects were negligible and included 
additional terms in the x- and z-momentum equations to account for centrifugal and Coriolis forces 
in the boundary layer due to blade rotation. These equations were coupled by using a modified 
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tangency boundary condition. This modified boundary condition enforces a velocity component 
(“transpiration velocity”) normal to the surface which “deflects the inviscid flow from the body 
surface thus simulating the displacement of the inviscid flow due to the momentum defect in the 
boundary-layer” [9]. Good comparisons were shown for integrated drag quantities (torque) on a 
non-lifting isolated rotor in hover for a range of hover tip Mach numbers. Also, good chordwise 
pressure coefficient comparisons were shown for two radial stations in a non-lifting, forward flight 
condition. Bridgeman, et ah [10] solved the unsteady, full potential equation coupled to a three 
dimensional boundary layer model for isolated rotors in hover and forward flight. This method 
is similar to that presented in [9], with a number modeling improvements. Though it is possible 
conceptually to include a fuselage in these full potential computations, this would be difficult in 
practice due to the blade-fixed coordinate systems typically used in such analyses. Thus, interac- 
tional aerodynamic computations would be difficult to compute using these existing tools. 

Recent solutions to the Euler equations for isolated rotor applications have used unstructured 
grid techniques [11, 12, 13] to refine the grid system efficiently to better capture wake structures 
such as tip vortices. These techniques may be extended easily to include a fuselage body. How- 
ever, since the Euler equations do not include viscous terms, computation of any viscous effects 
(e.g., viscous drag on a fuselage) would, like the full potential equation examples above, still re- 
quire coupling the method to a boundary layer analysis. 

Navier-Stokes computations have also been developed for use in rotorcraft analysis. While time 
averaged Navier-Stokes methods have been developed and are quite practical to use [14, 15, 16, 
17, 18], routine computations for unsteady flows on full configurations are not yet practical. Even 
so, several of these computations [19, 20] are found in the literature. Meakin [19] used a thin-layer 
Navier-Stokes method to calculate the flowfield around the Bell/Boeing V-22 Osprey tiltrotor air- 
craft, including the fuselage and the rotor, for a fictitious flight condition. Though this was a full 
aircraft simulation, the purpose of the computation was to demonstrate the feasibility of using a 
new domain connectivity algorithm for the moving, dynamic, overset grids for such a computa- 
tion. Since the computations were performed to demonstrate a technology, no comparisons are 
made with experimental quantities. Srinivasan and Ahmad [20] used a Navier-Stokes scheme to 
calculate the quasi-steady flowfield for a hovering rotor mounted on a whirl tower. Due to the 
quasi-steady nature of the hovering condition, the equations were solved in the blade-fixed coor- 
dinate frame, using a momentum source term in the equations to account for the centrifugal force 
of the blade rotation. This simulation was also a feasibility study, so only a comparison of the 
predicted and measured mean thrust values were presented. For this simulation, which utilized 
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approximately 1.3 million grid points, a computational time of 14 Cray-YMP hours was quoted. 
Ahmad and Duque [21] used a thin-layer Navier-Stokes method with embedded, moving, overset 
grids to demonstrate the ability to calculate the unsteady flowfield of an isolated, two-bladed rotor 
in forward flight. Blade surface pressures at several radial and azimuthal and local normal load co- 
efficients are compared to flight test data. These comparisons match reasonably well. Even though 
this computation did not include a fuselage body, the chimera grid scheme would render the task of 
including a body feasible, though at an additional computational cost. This isolated rotor compu- 
tation required substantial computational resources; it required approximately 45 Cray C-90 hours 
and generated 40 Gigabytes of data. 

Though most of these Navier-Stokes methods are suitable for computations over complete air- 
craft, including the helicopter rotor, all of the moving grid computations suffer from the require- 
ment to re-compute the grid domain connectivity information at every time step; this technique 
is known as a “dynamic chimera scheme” and has two distinct disadvantages. First, regenerating 
the grid domain connectivity at every time step can be as computationally expensive as, or even 
more computationally expensive than, the actual flow solution. In addition, in some instances, the 
time step is restricted not by the flow or flow solver, but by the moving grid domain connectivity 
requirement that a “hole point” not become a “field point” at any time step [19]. This requirement 
can potentially limit the time step not to physical phenomena, but to grid cell size. 

In addition to time step issues discussed above, other factors place limits on the current Navier- 
Stokes computations for rotorcraft. One issue is the numerical dissipation of concentrated vortices. 
It is well known that, in certain flight conditions, blade tip vortices can have a large effect on the 
rotor aerodynamics and that these vortices need to be computed in the flowfield over several rotor 
revolutions. Numerical studies discussed in the literature suggest that a 5th order scheme using 14 
points across the vortex core produces satisfactory results for a vortex that is well-aligned with the 
grid [22], However, in a typical rotorcraft simulation, the vortex location is not known a priori and 
thus the vortex in general will not be aligned well with the grid. Therefore a more strict resolution 
requirement is imposed on the numerical scheme [22]. For current methods either a prohibitively 
dense grid must be used to assure that the vortex is resolved well spatially, or grid adaption must be 
used to refine the grid in the regions that contain the vortices. Both methods are computationally 
expensive. 

Another such issue is turbulence modeling. Many of the turbulence models in current use were 
developed for wall bounded flows. They are not well suited to the three dimensional, non-isotropic 
turbulence associated with rotor blade tip vortices. 
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1.2.3 Hybrid Methods 

Considering the computational expense of current CFD methods, one possible approach to examin- 
ing unsteady rotor/fuselage interactional aerodynamics is to use hybrid methods. For unsteady ro- 
tor/fuselage aerodynamics, several hybrid methods have been developed. One such hybrid method, 
developed by Steinhoff, et al. [23], modified the Navier-Stokes equations by adding a “vorticity 
confinement” term to the momentum equations. This new term is used to prevent, or counter- 
act, numerical diffusion of concentrated vortical regions by “convecting” vorticity back toward the 
centroids of concentrated vorticity regions in the flowfield. This particular method is well suited 
for inclusion of a fuselage body. Boyd and Barnwell [24] first introduced a hybrid unsteady rotor 
model which weakly couples a Generalized Dynamic Wake Theory (GDWT) [25, 26, 27, 28] to a 
thin-layer Navier-Stokes model, OVERFLOW [29], Extensive induced inflow comparisons were 
made between Laser Doppler Velocimeter (LDV) measurements and predictions. Even though the 
computations were for an isolated rotor, excellent agreement was found with measured quanti- 
ties. Also presented was an outline of a method to couple a fuselage into the calculations using 
the overset grid capabilities in OVERFLOW. This new model uses the GDWT to obtain unsteady 
loading and unsteady induced inflow on the rotor, and then applies the unsteady loading inside 
OVERFLOW as a new unsteady pressure -jump, actuator disk-type, boundary condition. Another 
hybrid method, building on the previous literature [24], is developed in this research. 


1.3 Present Approach 

The objective of the current research is to develop an efficient, hybrid, unsteady computational 
model appropriate to the study of unsteady rotor/fuselage interactional aerodynamics. 

In examining fully CFD, unsteady, moving grid methods for complete rotorcraft, it can be ob- 
served that small time steps are needed for method stability, for capturing aerodynamic effects 
that are on the order of the rotor blade chord size, and for proper usage of the dynamic chimera 
grid scheme. However, to capture the primary effects of rotor/fuselage interactional aerodynamics, 
chordwise aerodynamics on the rotor blade are of less importance than the gross loading on the 
rotor blade itself. This can be seen by the successes of some of the singularity methods which use 
lifting line rotor blade models (Le., no chordwise loading distribution on the rotor blade) discussed 
previously in the “Literature Review” section above. In addition, fully CFD methods compute the 
rotor loading internally and require a number of rotor revolutions to obtain a periodic solution. The 
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combined requirements of needing very small time steps and of needing several rotor revolutions 
to obtain a periodic solution are a large contributor to the computational expense of these methods. 

A hybrid method is developed here which reduces the computational expense by separating 
the rotor loading calculation from the CFD component of the computation. This hybrid method is 
depicted schematically in figure 1.2. From the figure, it can be seen that there are three components 
to this hybrid method: 

1. Rotor Loading Model, 

2. Rotor/Fuselage Flowfield Model, 

3. Coupling Model. 

1.3.1 Rotor Loading Model 

The present approach separates the rotor loading model from the rotor/fuselage flowfield model. 
Splitting the procedure into these two separate models allows an otherwise computationally expen- 
sive element, the rotor loading computation, to be accomplished using a efficient, simplified model 
apart from the CFD computation. To compute the rotor loading, the GDWT [25, 26, 27, 28], which 
will be discussed in detail in a subsequent chapter, is used here. Previous implementations of the 
GDWT have focused on calculation of the unsteady inflow for an isolated rotor. As a significant 
advance over previous models, the unsteady rotor portion of the model uses the GDWT to calculate 
unsteady inflow and unsteady loading. The unsteady loading on the rotor is determined in the form 
of a A P, or “pressure jump”, across the rotor disk. This A P is then used as a boundary condition in 
the rotor/fuselage flowfield model. 

1.3.2 Rotor/Fuselage Flowfield Model 

The unsteady loading, as determined by the GDWT, is then used in conjunction with a Navier- 
Stokes model, in this case, OVERFLOW , to compute the time dependent flow over the fuselage 
including effects of a helicopter rotor. The loading is used in the Navier-Stokes method as an 
unsteady boundary condition in the flowfield. This boundary condition is effectively an unsteady 
actuator disk model. Though there will be concentrated regions of vorticity near the edges of 
an unsteady actuator disk model used as a boundary condition in this manner, these are not true 
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“tip vortices” and thus internal structure of these flow features is of secondary importance. This 
alleviates the need to develop prohibitively dense grids, use grid adaption, or use higher order 
schemes to resolve these concentrated vorticity regions. As such, turbulence modeling of the 
inside of these vortex structures becomes less important as well. So, by modeling the rotor as an 
unsteady actuator disk in OVERFLOW, several computationally expensive requirements typically 
needed for full CFD rotorcraft modeling are diminished. Using the above model, the Navier- 
Stokes method is then used to compute the periodic flowfield of the rotor/fuselage combination. 
Further details of this component will be discussed in a subsequent chapter. 

1.3.3 Coupling Model 

With the completion of the Navier-Stokes method, there are two solutions which were obtained 
with the same A P distribution: the GDWT solution, which is for an isolated rotor, and the OVER- 
FLOW solution, which contains both the rotor (as a boundary condition) and the fuselage body. 
The primary difference is that one solution contains a fuselage and the other does not. Since the 
loading in the GDWT depends on the rotor inflow, and since these inflow values are influenced by 
the presence of the fuselage, a method of coupling the two codes has been developed to account 
for the fuselage effects in the GDWT (and thus the rotor loading model). Since the fuselage effects 
on the rotor are assumed to consist of low frequency effects (as compared to the higher frequency 
blade-passage effects), the method employed here differences the time averaged inflow generated 
in the two successive solutions, and uses this difference as an “inflow correction” to the GDWT . 
This coupled process continues until only small differences are seen between successive iterations. 
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Computational Expense 


Figure 1.1: Analysis Types for Coupled Solutions. 
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Figure 1.2: Current Hybrid Method. 





Chapter 2 


Rotor Loading Model: Generalized 
Dynamic Wake Theory 

2.1 Introduction 

The Generalized Dynamic Wake Theory (GDWT ) of Peters, Boyd, and He [26], Peters and He 
[25], and He [27] is used in the present approach to obtain unsteady loading that is to be used in 
conjunction with OVERFLOW as discussed previously. In a later chapter, the OVERFLOW and 
the coupling technique between OVERFLOW and the GDWT will be discussed. Even though the 
GDWT is spelled out in detail in the literature, the present chapter will describe the GDWT for 
background purposes and describe the particular manner in which the theory is implemented for 
the current research. 


2.2 General Description 

The GDWT is a theory that was originally designed to pose the issue of unsteady aerodynamics of 
a helicopter rotor in a state-space form. This type of state-space form is desirable for inclusion in 
a rotor stability analysis since stability analyses for the rotor dynamics are usually presented in a 
state-space form as well. With the aerodynamics and dynamics of the rotor stated in similar forms, 
the solution of the system of equations is simplified. 
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2.3 GDWT Outline 

There are several aerodynamic concepts which are combined in the development of the GDWT. 
First, the acceleration potential derived by Kinner [30] for circular wing planforms is used with 
slight modification. The original acceleration potential derived by Kinner is a general form for the 
solution to Laplace’s equation (inviscid, linear, potential flow) in ellipsoidal coordinates. These 
modifications to Kinner’s acceleration potential function (or pressure function), applied by Peters, 
Boyd, and He [26], Peters and He [25], and He [27], are made to eliminate terms in the potential 
that are not compatible with boundary conditions associated with a rotor. This modified accel- 
eration potential, O, is then expressed using Legendre functions and transcendental functions in 
ellipsoidal coordinates as follows: 


<D(v,ti ,?,D = ” £ X ^r(v)Q7(iTl) K c (Dcos(m\j/)+ C(*) «"(«?)] t 2 - 1 ) 

^ m=0 n=m-b 1 

where the x™ and x™ terms are general coefficients of the pressure function and are, in general, 
functions of time and are determined from the loading on the rotor. The ellipsoidal coordinate 
system used here can be seen in figure 2. 1 . This figure shows a view of the xz plane with represen- 
tative r| and v values labeled. The $ coordinate (not shown in figure 2.1) is an angular, azimuthal 
coordinate, measured around the z-axis. The “rotor disk” is defined by the following conditions. 


( 2 . 2 ) 

(2.3) 

(2.4) 

where r is the radial coordinate on the rotor disk, measured from the axis of rotation, and \j/ is the 
angular, azimuthal coordinate measured about the axis of rotation. Equation (2.1) is effectively an 
expression for all admissible functions of loading on a rotor. To establish a link between loading 
on the rotor and induced inflow, the continuity equation and the linearized, incompressible Euler 
equations are used as follows: 


r\ — 0 
v = \f 1 — r 2 
? = V 




V 


\w' 


w 
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Si ~ v °°qi& = -®,i ( 2 - 6 ) 

qx = u (2.7) 

q 2 = v (2.8) 

<73 = w (2.9) 


( 2 . 10 ) 

where the summation convention is assumed over the index i, and £ is the coordinate pointing 
upstream along a streamline. Using equations (2.5) and (2.6), it can be shown that the pressure 
function, <I>, satisfies Laplace’s equation. Also, since these equations are linear, the pressure func- 
tion can be split into a component associated with each of the two terms on the left hand side of 
equation (2.6). Each of these components, in turn, also satisfies Laplace’s equation. As such, solu- 
tions can be derived for each component, then combined into a complete solution. The quantities 
used here are assumed to be total quantities, not perturbations. Also, it is assumed that only the 
velocity normal to the rotor disk is of interest and that it is of the following form: 

»’(?,¥.<) = £ £ $5W [<*5(0 co S (rH<)+P;(f)sin(^)] (2.11) 

r - 0 y==rf l,r+3,... 

With these assumptions, a closed form set of first order, non-linear differential equations, in state- 
space form, can be established for the induced inflow coefficients dc These equations are as 
follows: 


su 





(2.12) 


(2.13) 


With these equations, the unsteady aerodynamics and induced inflow are cast in the time domain 
and in a state-space form. The problem has now been reduced to the computation of the states of 
the model, a”, (3™, given a loading on the rotor. This is a non-linear model since the loading and the 
inflow are coupled through mass flow parameter, V™, in equations (2.12) and (2.13). As described 
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in reference [26], the mass flow parameter V™, which accounts for the energy added to the flow 
by the rotor, takes the place of the VU term that results from equation (2.6) in order to extent the 
theory and to cover the case of hover, where approaches zero. 


2.4 Solution Procedure 

Though the details of the GDWT are discussed in the literature, few details are provided for the 
solution of this set of first order, non-linear differential equations. The procedure used in the 
current research is described here. 

First, for a given set of rotor collective, lateral, and longitudinal pitch controls settings, the blade 
loading may be determined by any theory that can generate a loading given velocity (inflow) infor- 
mation. In the current research, as was done in previous literature [26, 25, 27], a two dimensional 
strip theory is employed. At first, this may seem to be an unnecessary restriction to a two dimen- 
sional theory. It has been shown [26, 25, 27] that this is not a restriction since the inflow and the 
loading are coupled. So, three dimensional effects, such as load reduction at the blade tip, are 
included in this theory. The equations of the strip theory used in this research are outlined below: 


0(r,\|/) = 0, w (r)+ 0 O + 0 C cos \i/ + 0 5 sin\jr 

(w + p sincts + Po^ cos V - O-5c0 + AX,) 

c t e ff = 0-- — : 

JJ r + p sinty 

L _ 7tc(r+ ^sin\|/) 2 a c // 
pQ 2 /? 3 ~ vT -M 1 

Equation (2.16) gives the lift at a particular point on the rotor disk. Through the effective angle 
of attack, a e ff, this loading theory includes local effects of the blade pitch (0), of the total induced 
inflow velocity (w), of the inflow due to the shaft tilt (^sinot*), of the inflow due to the mean blade 
coning angle (PoA* cost}/), of the velocity at the 3/4 chord point due to blade pitch rate (O.5c0), and 
of the inflow correction determined by the coupling scheme (AX;), described in a later chapter. In 
addition, the lift curve slope is assumed to be 2n per radian, the Prandtl-Glauert compressibility 
correction is applied to the lift, and a simple stall model is used which limits the maximum angle 
of attack to 10 degrees. Now, given the blade pitch settings of the rotor and the rotor operat- 
ing conditions, the local sectional loading can be determined from equations (2.14), (2.15), and 


(2.14) 

(2.15) 

(2.16) 



David D. Boyd, Jr. 


Chapter 2. Generalized Dynamic Wake Theory 


16 


(2.16). With the lift distribution determined, equations (2.12) and (2.13) are solved using a 4-stage 
Jameson-style Runge-Kutta technique [31] until periodic induced inflow is obtained; this usually 
occurs within two rotor revolutions. No blade dynamics model is used in this study and the blade 
hinge offset is assumed to be at the center of rotation ( i.e ., the flap hinge is at the center of rotation). 

With the above calculation complete, the mean thrust coefficient and mean hub moment coeffi- 
cients can be determined from the following: 


C T = 

4 / 
\/3n 

r-2n 

Cm — 

4 

3\/57t. 

rln 

/ x \ C d v 
Jo 

Cl = 

4 

3y/5n- 

rln 

/ \f 

Jo 


(2.17) 

(2.18) 

(2.19) 


Since typically the initial pitch settings do not produce the desired thrust and moment coef- 
ficients on the rotor, they must be adjusted in subsequent iterations until the desired thrust and 
moment coefficients are obtained. This is known as “trimming” the isolated rotor; this accounts 
only for the desired loading on the isolated rotor and does not include any “feedback” forces from 
any other source (such as a fuselage). The trim procedure used in the current research for the 
isolated rotor is a modified Newton-Raphson technique adapted from the literature [32], The fol- 
lowing equation is used to determine the pitch setting “corrections” which are used to iteratively 
trim the isolated rotor: 





( 2 . 20 ) 


where the matrix of partial derivatives is called the “derivative matrix” and each partial derivative 
is determined by using a one-sided forward difference formula and by using an independent per- 
turbation to each of the initial pitch settings. Each row in the derivative matrix is computed by 
an independent perturbation of each corresponding pitch setting and solving equations (2.12) and 
(2.13). Once computed, the derivative matrix is held unchanged throughout the subsequent rotor 
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trimming process. The ACy, ACa/, and A Cl are the changes in the current thrust and moment coef- 
ficients required to match the desired values. The trimming process is considered complete when 
value of the following is true: 


^(A Cr) 2 + (A Cm ) 2 + (AC t ) 2 < £ (2.21) 

where 6 is a specified tolerance. Now, with the trim task complete, the unsteady loading and 
unsteady induced inflow are known. 

At the end of the trim process, the unsteady loading is in the form of a sectional loading 
(i.e., force per unit span). However, for use in OVERFLOW, the loading needs to be in the form of 
a pressure that will be applied to a grid point in the flowfield. The sectional loading is converted 
to a pressure (force per unit area) using the assumption that the force is evenly distributed over 
the chordwise and spanwise extent of the local blade section of interest. These pressures are now 
ready for use in OVERFLOW, which will be discussed in a later chapter. 

Since the current research employs the GDWT in a manner in which it is not normally used, 
new computer coding has been developed here to compute required quantities from this theory, in 
combination with solution methods that have not previously been used with the GDWT. To validate 
that the new theory and models have been implemented correctly, a validation study is presented 
in the next chapter. 
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Figure 2.1: Ellipsoidal Coordinates. 









Chapter 3 


GDWT Validation 


3.1 Introduction 

Though there are numerous examples of the GDWT published in the literature, the present combi- 
nation of the GDWT with the present solution procedures has not been explored in the past. Also, 
computer coding to do the actual computations with these combinations of theory and solution pro- 
cedures was not available for the purposes of the current research. As such, validation is required 
to determine that the current model matches previously published literature on the subject. This 
validation effort is described in this chapter. 


3.2 Experimental Setup 

All of the experimental data used for the validation effort in this chapter was obtained from laser 
velocimeter inflow measurements made in the NASA Langley Research Center 14- by 22-Foot 
Subsonic Tunnel. The experimental setup and data used for this chapter is described in detail in 
references [33, 34, 35, 36]; for completeness, the experiments are described here in brief. 

Two different planforms were used in this test, one rectangular and one tapered. A summary of 
the two geometries is given in table 3.1 and a picture of the configuration, installed in the tunnel, is 
shown in figure 3.1. In this figure, the rotor/fuselage configuration consists of a generic helicopter 
fuselage body, known as the ROtor Body INteraction (ROBIN) fuselage [37], and one of the two 
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rotor configurations, described above. Inflow data was taken on these two rotor configurations at 
several advance ratios {i.e., forward speed divided by the rotor tip speed) for 12 azimuthal stations 
located 30 degrees apart (starting from an azimuth location of 0° , directly downstream of the rotor 
hub) and at a number of radial stations at each of these azimuth stations. Data samples were 
processed at a resolution of 128 samples per rotor revolution at each measurement location. For 
this chapter, the data is presented in two formats. The first format is the time averaged data, which 
is obtained by temporally averaging the unsteady data at each measurement location. The second 
format is time accurate data, which is the unsteady data before it is time averaged. 

There are two advance ratios used in this chapter. The first, p = 0.15, can be though of as 
boundary between (1) very low speed flight, where the fuselage flowfield is completely dominated 
by the rotor wake, and (2) moderate forward flight, where the fuselage flowfield is no longer 
dominated by rotor wake interactions. The second, p = 0.23 can be thought of as a moderate 
forward flight case. At both speeds, pressure pulse effects on the fuselage from the passing rotor 

blades is felt. 

The time averaged data will be presented in the form of contour plots of induced inflow ratio 
{i.e., magnitude of induced inflow velocity divided by rotor tip speed) mapped over the entire rotor 
disk. The time accurate induced inflow ratio will be presented as time histories which depict the 
induced inflow ratio at a particular fixed point over the rotor disk. For the contour plots of the 
measured data, it is worth noting that, even though the temporal resolution in the time accurate 
measurements corresponds to approximately 2.8° of blade travel, the spatial resolution in the mea- 
sured data is limited. This spatial resolution limit is imposed by the limited number of azimuthal 
locations at which measurements were taken over the rotor disk. For the contour plots, the impli- 
cation is that, since there were only twelve azimuthal measurement locations, the spatial frequency 
content of the induced inflow data in the azimuthal direction is limited by the Nyquist criterion to 
six harmonics per rotor revolution. In addition, there are some measurement locations at which no 

data is available. 

For the sole purpose of making consistent contour plots of measured data, which are presented 
subsequently, the time averaged data for induced inflow ratio over the rotor disk has been linearly 
interpolated onto a grid consisting of seventeen evenly spaced radial stations from the 20% ra- 
dius location to the tip location and onto twelve evenly spaced azimuth locations around the rotor 
disk. Interpolation onto a regular grid in this manner serves to “full-in” measured data that is not 
available at particular locations over the rotor disk. 
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Table 3.2 lists the matrix of cases, measured and predicted, that a will be used for the validation 
effort in this chapter. 


3.3 Time Averaged Induced Inflow 


The following two sections present the time averaged, measured and predicted induced inflow ratio 
for several configurations. The quantities presented here are ones that have been derived from the 
measurements and predictions by time averaging the time accurate induced inflow ratio. Also, 
it should be noted here that the measured quantities include the effects of the ROBIN fuselage 
geometry and the predicted quantities do not include these effects (i.e., the predictions are for an 
isolated rotor). 

3.3.1 Rectangular Planform 

Figure 3.2 shows a spatial contour plot comparison between the time averaged, measured and 
predicted induced inflow ratio for the rectangular planform at an advance ratio of 0.15, for a range 
of harmonics in the GDWT. Figure 3.2a is the measured data, including effects of the fuselage body. 
Figure 3.2b shows the GDWT predicted data using only two harmonics. Here, only the nominally 
linear streamwise gradient of induced inflow is predicted. Figure 3.2c shows the GDWT predicted 
data using four harmonics. In this figure, it can be seen that the major features of the measured 
data are predicted well. For example, the upwash on the forward section of the disk is predicted, 
though it does not cover as much of the forward region of the disk as in the measured data. The 
aft downwash regions, concentrated in the first and fourth rotor quadrants (see figure 3.18 for 
quadrant definitions), are also predicted. Figure 3. 2d shows the GDWT predicted data using eight 
harmonics. In this figure, all of the same features are present as in figure 3.2c, but the details 
of the induced inflow are slightly different. Figure 3.3 shows the same data as in figure 3.2, but 
these are lateral and longitudinal subsets of the data and show the radial variation of measured 
and predicted induced inflow. Figure 3.3a shows measured and predicted lateral variation of the 
induced inflow ratio for 2, 4, and 8 harmonics; the advancing and retreating sides of the rotor disk 
are labeled. Figure 3.3b shows measured and predicted longitudinal variation of the inflow ratio 
for 2, 4, and 8 harmonics; the forward and aft portions of the rotor disk are labeled. Figures 3.4 
and 3.5 show the same case as the previous two figures, except these predictions use 128 azimuth 
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steps per revolution instead of 64. No significant differences are apparent between the two different 
azimuthal resolutions; this shows that, for this particular case, 128 azimuth steps per revolution is 
a sufficient azimuthal resolution. 

Figures 3.6 to 3.9 show the same contour plots and lateral/longitudinal plots as in figures 3.2 to 
3.5, but for the higher advance ratio of 0.23. For this advance ratio, as with the lower advance ratio, 
little difference is seen between the use of 4 and 8 harmonics and little difference is seen between 
the use of 64 azimuthal steps and 128 azimuthal steps. 

3.3.2 Tapered Planform 

Figures 3.10 to 3.17 represent the same cases as in figures 3.2 to 3.9, but instead, using the tapered 
planform rotor. In comparison to the computations on the rectangular planform, the extent of 
the upwash region on the forward portion of the disk is better predicted based on the contour 
plots. Also, from the lateral and longitudinal plots, it can be seen that fewer harmonics are needed 
to predict the induced inflow behavior at the tip for the tapered planform. These findings for 
the contour plots and for the lateral and longitudinal plots are in agreement with findings in the 
published literature [25]. 


3.4 Time Accurate Induced Inflow 

The previous sections presented the time averaged quantities derived from the experimental and 
predicted data. The following two sections present the unsteady induced inflow ratio components 
associated with the cases presented above except for the cases corresponding to two harmonics. 
These two harmonic cases are excluded at this point because using two harmonics does not ade- 
quately represent the time averaged induced inflow ratio distribution over the rotor disk. In order 
to show the time dependent quantities more clearly, the local time averaged quantities have been 
removed from each of the following plots. Figure 3.18 shows the measurement locations where 
the comparisons of the unsteady induced inflow will be made. The radial and azimuthal positions 
are shown for locations A through J. Locations A through I are used in this chapter. Location J 
will not be used in this chapter. It is, however, used in chapters 7 and 8 and is shown here for later 
reference. 
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3.4.1 Rectangular Planform 

Figure 3.19 shows the unsteady component of induced inflow ratio at particular points on the rotor 
disk. The black dotted curves represent the measured data and the solid red lines represent the pre- 
dicted data. In many of the measured data, a distinct set of four pulses per revolution can be seen. 
These pulses represent blade passages past the measurement location. The predictions at these 
measurement locations also show a periodic signature that generally matches the measurements in 
phase of the signals. Comparing figure 3. 19 to figure 3.20 shows that for this configuration and this 
azimuthal resolution, the eight harmonic prediction matches the phase and amplitude better than 
the four harmonic prediction. Figures 3.21 and 3.22 show that the same holds for the 128 azimuth 
steps per revolution case. Figures 3.23 through 3.26 show that these results also hold for the higher 
advance ratio of 0.23. In these cases, it can be seen that, even using eight harmonics, the sharp, 
high frequency waveform of the measured inflow pulses is not matched well. 


3.4.2 Tapered Planform 

Figures 3.27 through 3.34 show the same features and results for the tapered blades as was pre- 
sented for the rectangular blades above. As with the rectangular blade results, there is a typical, 
high frequency four per revolution pulse (waveform) indicative of blade passages seen at the mea- 
surement locations; the phase of these pulses generally matches the phase of the measured data. 
Also, as with the rectangular blade predictions, the eight harmonic results predict the amplitude 
better than the four per revolution results. 


3.5 Observations 

From the evidence presented in this chapter, the following observations can be made. 

• Even though completely different methodology and computer coding from previous litera- 
ture was used in the implementation and solution of the equations of the GDWT, the current 
results match well with the previously published literature [25, 27]. 

• The best overall choice for number of harmonics, considering both the advance ratio and 
configuration variations, is eight. This conclusion is influenced more by the time averaged 
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predictions than the time accurate predictions. This “weighting” toward the time averaged 
computations is influenced by the fact that the coupling model (discussed later in chapter 5) 
uses the time averaged quantities, not the time accurate quantities. 

• Even with eight harmonics, it is not possible to simulate the sharpness of the blade passages. 

• The induced inflow ratio predicted by the GDWT is relatively insensitive to the number of 
azimuthal steps used. 

Based on the study presented in this chapter, and the observations presented in this section, it 
has been established that the new GDWT implementation is able to reproduce measured data in 
a manner comparable to previous literature. Thus, referring back to figure 1.2, the rotor loading 
model has been established. As such, the values of the A P pressure jump are now ready to be used 
in the rotor/fuselage model, which will be discussed in the next chapter. 
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Table 3.1: Rotor Geometries 


Property 

Rectangular 

Tapered 

radius 

0.8606 meters 

0.8255 meters 

root chord 

0.0660 meters 

0.0800 meters 

tip chord 

0.0660 meters 

0.0254 meters 

taper ratio 

(no taper) 

3:1 past 0.75R 

number of blades 

4 

4 

root cutout location 

0.24R 

0.24R 

flap/lag hinge location 

0.06R 

0.06R 

sweep of quarter-chord 

(none) 

(none) 

airfoil section 

NACA 0012 

NACA 0012 

twist 

-8° 

-13° 

precone 

(none) 

(none) 

nominal thrust coefficient 

0.0065 

0.0065 

solidity 

0.0977 

0.0977 

nominal tip Mach number 

0.55 

0.55 

approx, mean coning angle 

1° 

1° 

shaft tilt 

3° nose down 

3° nose down 
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Table 3.2 : Test and Prediction Matrix 


Advance Ratio 

No. Azimuths 

No. Harmonics 

Configuration 

0.15 

64 

2,4,8 

rectangular 

0.15 

128 

2,4,8 

rectangular 

0.23 

64 

2,4,8 1 

rectangular 

0.23 

128 

2,4,8 

rectangular 

0.15 

64 

2,4,8 

tapered 

0.15 

128 

2,4,8 

tapered 

0.23 

64 

2,4,8 

tapered 

0.23 

128 

2,4,8 

tapered 



David D. Boyd, Jr. 


Chapter 3. GDWT Validation 


27 



Figure 3.1: ROBIN Fuselage in the NASA Langley Research Center 14- by 22-Foot Subsonic 
Tunnel(1986). 


David D. Boyd, Jr. 


Chapter 3. GDWT Validation 


28 


measured 



downwash 

upwash 

predicted 
2 harmonics 



predicted 
4 harmonics 



predicted 
8 harmonics 



Figure 3.2: Measured and GDWT Predicted Time Averaged Induced Inflow Ratio; Rectangular 
Planform, u = 0. 15, 64 Azimuths. 
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Figure 3.3: Measured and GDWT Predicted Lateral and Longitudinal Time Averaged Induced 
Inflow Ratio; Rectangular Planform, p = 0. 15, 64 Azimuths. 
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Figure 3.4: Measured and GDWT Predicted Time Averaged Induced Inflow Ratio; Rectangular 
Planform, u. = 0. 15, 128 Azimuths. 
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Figure 3.5: Measured and GDWT Predicted Lateral and Longitudinal Time Averaged Induced 
Inflow Ratio; Rectangular Planform, p = 0.15, 128 Azimuths. 
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Figure 3.6: Measured and GDWT Predicted Time Averaged Induced Inflow Ratio; Rectangular 
Planform, u = 0.23, 64 Azimuths. 
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Figure 3.7: Measured and GDWT Predicted Lateral and Longitudinal Time Averaged Induced 
Inflow Ratio; Rectangular Planform, p = 0.23, 64 Azimuths. 
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Figure 3.8: Measured and GDWT Predicted Time Averaged Induced Inflow Ratio; Rectangular 
Planform, p = 0.23, 128 Azimuths. 
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Figure 3.9: Measured and GDWT Predicted Lateral and Longitudinal Time Averaged Induced 
Inflow Ratio; Rectangular Planform, p = 0.23, 128 Azimuths. 
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Figure 3.10: Measured and GDWT Predicted Time Averaged Induced Inflow Ratio; Tapered Plan 
form, ju = 0. 15, 64 Azimuths. 
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Figure 3.11: Measured and GDWT Predicted Lateral and Longitudinal Time Averaged Induced 
Inflow Ratio; Tapered Planform, n = 0. 15, 64 Azimuths. 
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Figure 3.12: Measured and GDWT Predicted Time Averaged Induced Inflow Ratio; Tapered Plan 
form, u = 0.15, 128 Azimuths. 
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Figure 3.13: Measured and GDWT Predicted Lateral and Longitudinal Time Averaged Induced 
Inflow Ratio; Tapered Planform, p = 0.15, 128 Azimuths. 
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Figure 3.15: Measured and GDWT Predicted Lateral and Longitudinal Time Averaged Induced 
Inflow Ratio; Tapered Planform, /u = 0.23, 64 Azimuths. 
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Figure 3.16: Measured and GDWT Predicted Time Averaged Induced Inflow Ratio; Tapered Plan- 
form, p = 0.23, 128 Azimuths. 



David D. Boyd, Jr. 


Chapter 3. GDWT Validation 


43 


predicted, 2 harmonics 

predicted, 4 harmonics 

predicted, 8 harmonics 

o o o o o measured 



r/R (a) r/R (b) 


Figure 3.17: Measured and GDWT Predicted Lateral and Longitudinal Time Averaged Induced 
Inflow Ratio; Tapered Planform, p = 0.23, 128 Azimuths. 
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Figure 3.19: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Rectangular Planform, p = 0.15, 64 Azimuths, 4 Harmonics. 
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Figure 3.20: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Rectangular Planform, p = 0. 15, 64 Azimuths, 8 Harmonics. 
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Figure 3.21: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Rectangular Planform, p = 0.15, 128 Azimuths, 4 Harmonics. 
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Figure 3.22: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Rectangular Planform, p = 0. 15, 128 Azimuths, 8 Harmonics. 
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Figure 3.23: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Rectangular Planform, p = 0.23, 64 Azimuths, 4 Harmonics. 
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Figure 3.24: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Rectangular Planform, p — 0.23, 64 Azimuths, 8 Harmonics. 
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Figure 3.25: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Rectangular Planform, p — 0.23, 128 Azimuths, 4 Harmonics. 
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Figure 3.26: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Rectangular Planform, p = 0.23, 128 Azimuths, 8 Harmonics. 
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Figure 3.27: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Tapered Planform, p = 0.15, 64 Azimuths, 4 Harmonics. 
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Figure 3.28: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Tapered Planform, p = 0. 15, 64 Azimuths, 8 Harmonics. 
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Figure 3.29: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Tapered Planform, p = 0.15, 128 Azimuths, 4 Harmonics. 
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Figure 3.30: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Tapered Planform, p = 0.15, 128 Azimuths, 8 Harmonics. 


Inflow Ratio - Mean Inflow Ratio • Mean Inflow Ratio - Mean 


David D. Boyd, Jr. 


Chapter 3. GDWT Validation 


57 


measured 

predicted 



'F = 30° 



oo5 r Location C 




-0 05 »-*- *»*l«i.**it**iilt t ■* * I 

0 90 180 270 300 


¥=180° 




'P = 300° 



Ref. Blade Location 


Ref. Blade Location Ref. Blade Location 


Figure 3.3 1 : Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Tapered Planform, jj = 0.23, 64 Azimuths, 4 Harmonics. 
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Figure 3.32: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Tapered Planform, p = 0.23, 64 Azimuths, 8 Harmonics. 
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Figure 3.33: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Tapered Planform, ju = 0.23, 128 Azimuths, 4 Harmonics. 
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Figure 3.34: Measured and GDWT Predicted Unsteady Induced Inflow Ratio With Mean Values 
Removed; Tapered Planform, p - 0.23, 128 Azimuths, 8 Harmonics. 





Chapter 4 

Rotor/Fuselage Flowfield Model: 
OVERFLOW 
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4.1 Introduction 


A rotor loading model which computes a distributed A P pressure jump value on the isolated rotor 
disk has been established and discussed extensively in chapters 2 and 3. From figure 1 .2, it can be 
seen that these A P values are used in the rotor/fuselage flowfield model to compute the entire flow- 
field of the combined rotor/fuselage system. The rotor/fuselage flowfield is computed by solving 
the Navier-Stokes equations in conjunction with conservation of mass and energy. The Reynolds 
averaged Navier-Stokes equations in Cartesian coordinates are as follows: 


where. 


dQ gcg-fl) 9(3- Gy) a (ff-fly) 

dt dx dy dz 
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In the above set of equations, p is the density, u,v,w are the velocity components, p is the 
pressure, eo is the total energy per unit volume, T j are components of the stress tensor, and q x ,qy, q z 
are components of the heat flux. 

The quantities above are all considered to be averaged, or mean-flow, quantities. Also, Newto- 
nian flow is assumed. Thus, the stress tensor is proportional to the velocity gradients in the flow. 
The constant of proportionality is composed of the sum of two components. The first component 
is the molecular viscosity, pi, which accounts for laminar or molecular viscosity, and is a property 
of the fluid. The second component is the turbulent eddy viscosity, p t , which accounts for the 
turbulence in the flow and is computed using a turbulence model. 

The above system of equations is then closed by introducing the perfect gas law as follows: 


P = P%T 


(4.6) 


where is the gas constant and T is the temperature. 
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For the current rotor/flowfield model, a time accurate, Reynolds averaged Navier-Stokes (RANS) 
tool known as OVERFLOW [29] is used to solve the above set of equations. OVERFLOW is a 
computer code that uses a finite difference technique to solve the compressible, RANS equations 
in generalized coordinates. Even though OVERFLOW is technically a RANS solver, it is typically 
used in a thin-layer mode by ignoring the viscous terms that are not associated with the direc- 
tion normal to surfaces in the flowfield; that approach is used here as well. The code employs 
a chimera (overset) grid scheme [38], which is helpful for complex geometries. Several solution 
procedure options are available in this code. A second or fourth order central difference scheme 
with a second and fourth order artificial dissipation scheme is used for the convective and viscous 
terms. A Pulliam-Chaussee scalar diagonal inversion [39] is typically used for the left-hand of the 
equations, though other options are possible, namely a Beam- Warming [40] scheme and a Lower 
Upper Symmetric Gauss Seidel (LU-SGS) scheme [41]. For steady state computations, a local 
time stepping scheme and a multigrid scheme [42] are implemented to accelerate convergence. 
Also, for steady state computations at low Mach numbers, a low Mach number preconditioning 
scheme [42] is used. For unsteady computations, a “Newton sub-iteration” scheme [43] is im- 
plemented for each time step to reduce linearization and factorization errors and to increase the 
time accuracy of the scheme to approximately second order. Also available are several turbulence 
models and an extensive set of boundary condition options. Available turbulence models include a 
Baldwin-Lomax model, a Baldwin-Barth model, a Spalart-AIlmaras model, and a k — co model; the 
Spalart-Allmaras turbulence model is used for all computations presented here. Boundary condi- 
tion options include conditions for inviscid walls, viscous walls, periodic grids, symmetry planes, 
singular axes, inflow, outflow, characteristic conditions, etc. For both steady and unsteady com- 
putations, all boundary conditions are applied explicitly; some of these boundary conditions may 
also be applied to any region inside the volume of the computational grid, not just at the outer grid 
boundary faces. 

This chapter describes the manner in which OVERFLOW is used for the rotor/fuselage flowfield 
model depicted in figure 1.2. 

4.2 New Boundary Conditions 

For the purposes of the current research, OVERFLOW has been extended to include two new, 
novel, explicit boundary conditions which use the pressure jump previously computed by the 
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GDWT to model a helicopter rotor. One of the boundary conditions is used for time averaged 
computations and the other is used for unsteady, time accurate computations. Both are applied to 
two planes of a non-rotating, cylindrical grid with an “iblanked plane” in between. Figure 4.1a 
shows a top view schematic of the rotor disk used in the new boundary conditions. This schematic 
represents the non-rotating, cylindrical portion of the grid that is used to represent the rotor disk. 
This “rotor portion” of the grid is a subset of the much larger overset volume grid set that is used to 
represent the entire rotor/fuselage flowfield. This grid subset is where the new boundary conditions 
are applied. 

The shaded area in this figure represents a small region of the cylindrical grid; the rectangle 
outlined by the dark lines represents a rotor blade. Both of these regions are enlarged in the figure 
for clarity of presentation. Figure 4.1b shows an edge view of this same schematic. In this view, 
five planes can be seen. These planes consist of an upper rotor plane, an iblanked plane, a lower 
rotor plane, and two planes, labeled A and B, which are used in the formulation of the boundary 
conditions. 

Although the focus of the current research is on unsteady interactions, it is necessary to have 
both the time averaged and time accurate boundary conditions discussed above. These boundary 
conditions are used in a complementary manner as follows. First, the steady state flowfield around 
the isolated fuselage is determined. Then, using this isolated fuselage computation as a starting 
point, the steady state flowfield of the fuselage, including the time averaged rotor model, can be 
determined using the time averaged rotor boundary condition. The previous two stages are used 
to set up the steady state flow features so that the unsteady computations can be used as a final 
stage, which is executed until a periodic solution is obtained. This “building block” approach of 
using a steady state computation as an initial condition to a time averaged computation, and in 
turn, using the time averaged computation as an initial condition to the time accurate computation, 
is used to reduce the computational resources required relative to using a completely time accurate 
computation from inception. 

Next, the details of each of the new boundary conditions will be described. 
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4.2.1 Time Averaged Boundary Condition 

As stated above, the boundary conditions are applied on two planes, separated by an “iblank 
plane 1 ”. Use of the “iblank plane” between the upper and lower rotor planes prevents the artif- 
ical dissipation in OVERFLOW from being activated by the new pressure jump imposed across 
the two planes. 

The steps in the application of this boundary condition are as follows: 

1. Identify the two planes (labeled A and B in figure 4.1b) surrounding the upper and lower 
rotor planes, 

2. Average the conservative flow quantities in planes A and B as follows: 

Qi„ = (4.7) 

where Qi are given in equation (4.2), 

3. Replace the existing conservative flow quantities in the upper and in the lower rotor planes 
with these average values, 

4. Add an “additional energy term” to the fifth conservative flow quantity (Qs). 

Steps 1, 2, and 3 above are relatively straightforward and are accomplished with application 
of existing boundary conditions. Step 4 is the new step in this process introduced here and is 
described below. 

As stated above, both boundary conditions use the pressures from the GDWT with a few slight 
modifications. The first new boundary condition time averages the unsteady pressures at each radial 
and azimuthal location, multiplies by the ratio of the local blade area to the local computational cell 
area, A(r), to maintain the same level of thrust between the two methods, and divides by (y— 1) to 
convert the pressures into an energy-like term that is compatible with the Qs variable above. This 
conversion to a Q$ compatible quantity can be combined into one expression as follows: 

4(f) Nt 

[(P *0 )odd] (r, v) = v, t) (4.8) 

'The term “iblank” refers to a technique employed in many methods that use the chimera scheme. This technique 
involves intentionally excluding certain points (here, a certain plane) from the computation of the flowfield. 


David D. Boyd, Jr. 


Chapter 4. OVERFLOW 


66 


where Nt is the number of time steps used in one revolution. The value of the left side of equa- 
tion (4.8) is then split in half. One half of the term is added to the peo equation for the lower 
plane of the rotor grid at the current azimuth and radial location, the other half is subtracted from 
the p<?o equation for the upper plane in the rotor grid at the current azimuth and radial location; 
this effectively adds energy to the flowfield, to mimic an actuator disk. The splitting of the addi- 
tional peo term and placement on either side of an “iblank plane” prevents the artificial dissipation 
in OVERFLOW from acting on the effective additional pressure jump; the artifical dissipation in 
OVERFLOW is designed so that it will not operate across an iblank region. Without this iblank 
plane, the artificial dissipation would operate on the pressure jump, smoothing the pressure jump 
unnecessarily. 

4.2.2 Time Accurate Boundary Condition 

The new time accurate boundary condition is applied in a manner similar to that used for the 
time averaged boundary condition above. The differences are that (1) the pressures are not time 
averaged before they are converted to energy-like terms, and (2) the pressures are evenly distributed 
over all azimuthal grid lines that cross the blade chord at the local blade radial station. Similar to 
the time averaged boundary condition (with exception (1) above), the following additional energy- 
like term is determined: 


[(p*0 )odd](r, V, 0 = AP(r, y, t) (4.9) 

As in the time averaged boundary condition above, the value of the left side of equation (4.9) is 
then split in half. One half of the term is added to the peo equation for the lower plane of the 
rotor grid, the other half is subtracted from the peo equation for the upper plane in the rotor grid. 
As discussed above, splitting the additional term and applying it on either side of an iblank plane 
prevents the artifical dissipation term from acting on the pressure jump. Also, at each radial station 
on the blade, the number of grid lines that lie on the chord are computed based on the chord length 
and the arclength between successive azimuthal grid lines. The additional energy term above then 
is distributed evenly over these grid lines. This distribution allows the shape of the blade to be 
better modeled in the cylindrical grid. 


David D. Boyd, Jr. 


Chapter 4. OVERFLOW 


67 


4.3 Summary 

Referring to figure 1.2, a rotor/fuselage flowfield model has been introduced in this chapter. This 
new model consists of a Navier-Stokes method in which several new boundary conditions have 
been developed. These boundary conditions use information from the GDWT to model an addition 
of energy to the flowfield as in an actuator disk method, but in an time accurate setting. In the next 
chapter, the third box in figure 1.2, “coupling model”, will be discussed. 


David D. Boyd, Jr. 


Chapter 4. OVERFLOW 


68 


<w' 




A 





lower rotor plane 
"iblanked" plane 
upper rotor plane 


(b) 


Figure 4. 1 : Rotor Schematic for New Boundary Conditions. 




Chapter 5 


Coupling Model 


5.1 Introduction 

Since the pressure jump determined by the isolated rotor loading model depends on the inflow 
ratio distribution over the rotor disk, and since the presence of a fuselage alters this inflow ratio 
distribution, some type of coupling model is needed to adjust the inflow ratio and pressure jump 
distributions to reflect the effects of the fuselage. The current coupling method adopts an “inflow 
correction” model. In this type of model, the loading on the rotor is determined based on the inflow 
at the rotor disk, including additional inflow, or inflow corrections, generated by the presence of the 
fuselage. Presently, the forces and moments on the fuselage are not accounted for directly in the 
computation of the rotor forces; the rotor forces are, however, computed using inflow corrections 
to account for the presence of the fuselage. Since typically, for a model in a wind tunnel, the rotor 
is trimmed to a specified thrust, independent of the forces acting on the fuselage, this modeling 
assumption is valid here. 

In previous chapters, discussions of the rotor loading model and the rotor/fuselage flowfield 
model were presented. In this chapter, the lower box of figure 1.2, the “coupling model will be 
discussed. 
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5.2 Isolated Rotor Configurations 

Previous literature [24] used a hybrid method similar to that presented here. The current research 
builds on the idea presented in reference [24], which is a subset of the current work. This subset is 
contained inside the dashed box region in figure 5.1. For an isolated rotor computation, no coupling 
model, or feedback loop, is needed to provide fuselage correction information to the rotor loading 
model. Therefore, an isolated rotor computation is determined by a single pass through the method, 
stopping at the end of the rotor/fuselage flowfield computation. This was the method employed in 
reference [24]. 


5.3 Rotor/Fuselage Configurations 

When there is a fuselage present in the rotor/fuselage flowfield model, a correction scheme is 
needed to feed information back into the isolated rotor loading model to account for the fuselage 
presence. The direction of the arrows in figure 5.1 indicates that the correction scheme should 
reflect the effects of the fuselage on the rotor system, as opposed to the rotor effects on the fuselage. 

In observing effects of a fuselage on the rotor system, it can be noted that time scales associated 
with the flow around the fuselage are much larger than those associated with the rotation of the 
rotor system. The effect of the fuselage on the surrounding fluid is one of slow displacement. 
However, due to its rotation, the rotor blade will pass this fluid particle at a much higher speed. 
Thus, for a given time increment, the rotor blade will traverse far more distance than a fluid particle 
traveling along the fuselage. 

Based on these scale differences, it is expected that the presence of the fuselage produces effects 
on the rotor disk that are not spatially or temporally concentrated when viewed in the fuselage fixed, 
spatial frame of reference. However, when this effect is viewed in the rotor blade rotating frame 
of reference, the rotor blade “sees” the spatial distribution of inflow generated by the presence of 
the fuselage as a slowly varying (temporal) inflow quantity. This leads to the conclusion that the 
influence of the fuselage on the rotor system can be viewed as a time averaged perturbation 1 to the 
isolated rotor configuration in the fuselage fixed reference frame. 

With the view that the fuselage effect on the rotor system can be considered to be a time averaged 


'Note that here, a perturbation is not necessarily small. 
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inflow ratio correction distributed over the rotor disk, an efficient coupling method which accounts 
for the primary fuselage effects on the rotor inflow ratio, can be developed. The coupling method 
developed here is of that type. 

To determine the inflow corrections to be fed back into the rotor loading model, a difference is 
taken between the time averaged, filtered inflow ratio from the rotor/fuselage flowfield model and 
the time averaged inflow ratio from the GDWT: 

AX; = ?(lf FFM ) - (5.1) 

where jF(-) is a filtering operation performed on \f FFM . To show the reason for the filtering oper- 
ation on the Xf FFM term, it is necessary to examine the frequency limitations of each component in 
the method. For the rotor loading model component, the spatial frequency content of the induced 
inflow ratio determined by the GDWT, X FLM , is limited to the number of harmonics chosen in the 
method. For example, it was shown in chapter 3 that eight harmonics of rotor inflow are sufficient 
for that model in this context. Therefore inflow ratio information, time averaged or time accurate, 
is limited to eight harmonics (in this example). Note that this limitation does not imply that the 
loading distribution computed by the GDWT is limited to the given number of harmonics. This is 
because the GDWT computes the inflow ratio (up to a specified number of harmonics) given any 
loading distribution. This computed inflow ratio influences the loading distribution through the w 
term in equation (2.15). Thus, the otherwise two dimensional loading distribution reflects inflow 
corrections up to the number of harmonics in the model. 

While the GDWT limits the number of inflow harmonics computed to a specified value, the 
inflow computed by the rotor/fuselage flowfield model, Xf FFM , is only limited by the spatial reso- 
lution found in the rotor grid. For example, a typical rotor grid used in this method [24] uses 128 
azimuthal grid lines per rotor revolution. Using the Nyquist cutoff concept, this means that the 
limit of the frequency content in the azimuthal direction is 64 harmonics. If AX, were computed 
and used in the rotor loading model without the filtering operation, inflow corrections (and there- 
fore loading distribution corrections) would be made inside the GDWT that are inconsistent with 
the inflow corrections made internal to the model using the w term in equation (2. 15). To eliminate 
this inconsistency, the Xf FFM term is filtered to match the frequency content of the X FLM term. 

A consistent filtering operation, ?(•), is derived from the GDWT method and is given in Ap- 
pendix A. Applying this filtering, J : (%f FFM ), then computing AX, from equation (5.1) above, this 
term may be included in equation (2. 15) to complete the coupling loop. With this coupling method, 
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all of the boxes of figure 1.2 have been described. The next few chapters will discuss applications 
of entire method. 
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Figure 5.1: Current Hybrid Method. 





Chapter 6 


Results: Isolated Fuselage 


6.1 Introduction 

Chapters 1, 2, 4, and 5, have described the components of the new computational model for ro- 
tor/fuselage unsteady interactional aerodynamics and how they are used. In the next three chapters, 
a full configuration will be constructed from basic components: (1) an isolated fuselage, (2) an iso- 
lated rotor, and (3) a rotor/fuselage combination. The subject of this chapter is the isolated fuselage 
component. Even though, strictly speaking, an isolated fuselage model is not a new computational 
model, it is necessary to establish that the computations can be compared favorably to measured 
data. 


6.2 Experimental Setup 

The experimental setup for the comparisons in this isolated fuselage chapter are described in detail 
in Freeman and Mineck [37]. However, for completeness, the setup will be discussed here in brief. 

A test was conducted in the NASA Langley Research Center 14- by 22-Foot Subsonic Tunnel, 
in which steady state surface pressures were obtained at one hundred and sixty-two locations on 
a generic helicopter fuselage shape for a number of flight conditions with and without a rotor in- 
stalled. This fuselage shape, which is derived mathematically using “super-ellipse” equations [37], 
is known as the ROtor Body INteraction (ROBIN) fuselage (see figure 6.1). In the photograph 


David D. Boyd, Jr. 


Chapter 6. Results: Isolated Fuselage 


75 


shown in figure 6. 1, the ROBIN is sting mounted, with the rotor installed. For the isolated fuselage 
portion of the test, the blades were removed, but the hub was left in place. Though test condi- 
tions included configurations with and without the rotor installed, only steady state pressures were 
measured on the fuselage. As such, only the isolated fuselage pressures will be used here. 


6.3 Computational Grid System 

A system of thirteen grids, combined using the chimera grid capabilities in OVERFLOW, is used to 
represent the ROBIN fuselage geometry and flowfield. The sting mount, rotor hub, and wind tunnel 
walls are not modeled. A listing of all of grid names and sizes are displayed in table 6.1. From 
table 6.1 it can be seen that five of these grids involve surface grids associated with the fuselage. 
These surface grids are shown and labeled in figure 6.2. Note that, for clarity of presentation, only 
every third grid line is plotted in figure 6.2. All remaining grids are field volume grids. Figure 6.3 
shows a view of the grids immediately surrounding the fuselage. For clarity of presentation, every 
second grid line in this figure has been removed. Also seen is the positioning of the rotor grid 
with respect to the fuselage. For the isolated fuselage computations presented in this chapter, the 
rotor grid is simply another volume grid (i.e., there is no rotor boundary condition applied on the 
rotor grid). The far field grid extents approximately 2.5 fuselage lengths upstream of the fuselage 
nose, downstream of the fuselage tail, to the left and right of the fuselage, and above and below the 
fuselage. 

These grids were developed such that, with minimal effort and minimal grid modification, all 
configurations used in the current research (including an isolated fuselage, an isolated rotor, and 
a rotor/fuselage combination) could be modeled simply by "grid replacement ’. As an example, 
when converting from a rotor/fuselage configuration to an isolated rotor configuration, only the 
grids containing the fuselage need to be replaced; these are replaced by regular volume grids that 
maintain the outer boundary shape of the fuselage grids. All other grids remain the same. 

These grids were also designed to maintain a “double fringe” overlap between all overlapping 
grids (except for the outer field grid) to maintain the second and fourth order artificial dissipation 
scheme used in OVERFLOW. In addition, grid spacings were refined to maintain viscous spacings 
at the fuselage surface such that there is at least one grid point in the laminar sub-layer region of 
the boundary layer and to maintain approximately 1 .5 million total grid points in order to keep the 
computational expense reasonable. 


David D. Boyd, Jr. 


Chapter 6. Results: Isolated Fuselage 


76 


One novel feature of this grid system is the manner in which the rotor grid is included. In 
previous isolated fuselage studies of the ROBIN fuselage [44, 45] the grid system did not include 
a rotor grid, thus a hyperbolic grid generator was used to “grow” the volume grids out from the 
surface grids without regard to the location of the outer boundaries of these grids. However, due to 
the geometry of the outer boundaries of grids produced by a hyperbolic grid generator, inclusion of 
a cylindrical rotor grid by the technique of “hole cutting” [46], while maintaining a double fringe 
chimera scheme, is difficult. This problem arises due to the lack of control of the outer boundary 
shape of grids generated by hyperbolic grid generation methods. To alleviate this problem, an 
elliptic grid generation method is used to control the outer boundary shape of some grids. With 
this technique, the overlap of the grids can be better controlled to produce a double fringe overlap 
and the amount of “hole cutting” can be substantially reduced. 

For the current research, a hybrid set of grids generated by hyperbolic and elliptic grid generators 
[47 5 48, 49, 50] is used. For volume grids in which no control was needed over the outer boundary 
shape, a hyperbolic grid generator was used. If control over the outer boundary geometry was 
desired, an elliptic grid generator was used. 


6.4 Steady State Pressure Prediction 

As discussed previously, the first step in the OVERFLOW computations is to execute the isolated 
fuselage configuration in a steady state mode. For the current computations, the default settings 
recommended for OVERFLOW [29] are used, except where noted. These defaults include central 
differencing of the right hand side of the equations, scalar diagonal inversion of the left hand side 
of the equations, a matrix dissipation scheme, low Mach number preconditioning, multigrid, full 
multigrid (mesh sequencing), local time stepping with a minimum CFL number constraint, and the 
Spalart-Allmaras turbulence model. The case presented here is for an angle of attack of 0 and a 
freestream Mach number of 0. 1265. 

Figure 6.4 shows the force coefficients in the lift, drag, and sideward directions as a function 
of iteration number in OVERFLOW. These force and moment coefficients are defined as follows 
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where a list and derivation of reference quantities for these is given in reference [51]. Though these 
do not directly correspond to the standard lift, drag, and side force coefficients, they do provide 
a useful insight into the convergence behavior of the OVERFLOW iteration from an engineering 
standpoint. 

Figure 6.5 shows the pressure tap locations on the ROBIN fuselage. The locations are shown 
with respect to the downstream coordinate (^-direction). Along the top of the figure, the lowercase 
letters indicate the locations where predictions are compared to measurements in the subsequent 
figure. 

Figure 6.6 compares measured and predicted pressure coefficients vs the vertical location on 
the fuselage at various places stations along the length of the fuselage. The vertical axis in figure 
6.6 is plotted in the standard negative C p format and the horizontal axis is the vertical location 
(z coordinate) of the current section; the solid lines are the predictions and the symbols are the 
measured values. 

The pressure coefficient used here is defined in equation (6.4). From a practical standpoint, 
computation of C p inside OVERFLOW is accomplished using equation (6.5) and equation (4.2), 
taking into account the nondimensionalizations used in OVERFLOW. 


C p = 


P-P~ 

IPV* 


(6.4) 


Cp = 


2(7- 1) 
Ml 


Q 5 


@2 + @3 + @3 1 


201 


7(7-1). 
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In figure 6.6 it can be seen that predicted pressure coefficients match the experimental data well 
in most areas. At locations (f) and (g), discrepancies can be seen on the upper and lower portions 
of the fuselage. These differences could be due to the presence of the hub and sting mount in the 
experimental setup. In this figure, predicted C p values are plotted from both sides of the fuselage 
for this symmetric flight condition. In this figure, predicted pressure coefficients from both sides 
of the fuselage are plotted since the computation included the full fuselage instead of a half body 
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fuselage plus a symmetry condition. As, expected, no difference is seen between the two sides. 
In addition to matching well the experimental data, these predictions are consistent with similar 
predictions in the literature [45], which used other Navier-Stokes methods (codes), for the same 
configuration and flight condition. 


6.5 Observations 

OVERFLOW has the capability to predict the steady state pressure coefficients on an isolated 
fuselage configuration in an incompressible flight condition, and these predictions are consistent 
with similar computations in the literature. 
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Table 6.1: Computational Grid System 


Grid Name/Description 

Dimensions 

Fuselage Grid 

93 x 117x25 

Nose Grid 

30 x 117x33 

Tail Grid 

28 x 117x33 

Collar Grid 

29 x 1 1 1 x 25 

Top Grid 

41 x 41 x 25 

Lower Field Grid 

55 x 50 x 12 

Upper Field Grid 

55 x 50 x 22 

Nose Field Grid 

42 x 70 x 43 

Tail Field Grid 

30 x 39 x 29 

Left Field Grid 

55 x 32 x 30 

Right Field Grid 

55 x 32 x 30 

Outer Field Grid 

93 x 82 x 26 

Rotor Grid 

37 x 129 x 43 
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Figure 6. 1 : ROBIN Fuselage in the NASA Langley Research Center 14- by 22-Foot Subsonic 
Tunnel. 
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Figure 6.5: Pressure Tap Locations 
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Figure 6.6: Measured and Predicted Pressure Coefficients vs Vertical Location for a — 
Moo = 0.1265. 
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Chapter 7 

Results: Isolated Rotor 


7.1 Introduction 

Chapter 6 showed that the pressure distribution could be predicted on the isolated ROBIN fuselage 
configuration for a representative forward flight condition of a rotorcraft. This chapter will employ 
the isolated rotor method as described in chapter 5. Previous predictions using the current model 
[24] compared well with time averaged and time accurate laser velocimeter (LV) data. In that work, 
the predictions were for an isolated rotor, whereas the experimental data contained the effects of a 
fuselage. Subsequent to those predictions, new experimental data has been acquired on an isolated 
rotor configuration (described below). Comparisons between these data and the current method 
are the subject of this chapter. 


7.2 Experimental Setup 

For comparisons in this section, results from an isolated rotor configuration in the NASA Langley 
Research Center 14- by 22-Foot Subsonic Tunnel are shown. The isolated rotor test system (IRTS) 
[52, 53] is shown in figure 7.1. Here, it can be seen that the rotor is suspended from the ceiling 
of the tunnel from a tapered, cylindrical shaped sting. A three component laser velocimeter (LV) 
system was used to measure the three components of velocity at a point at an azimuth of 84 , a 
radial station of r/R = 0.81, and one blade chord above the tip path plane of the rotor. The rotor 
was trimmed to a nominal 0° flapping angle (e.g., the first flapping harmonics are Pi c = = 0 ), 
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a nominal thrust coefficient of 0.0064, a nominal shaft tilt, a s , of 3° nose down, and a freestream 
Mach number of 0.1265. The LV data was processed at 128 samples per rotor revolution. 


7.3 Computational Grid System 

As discussed in chapter 6 for the isolated fuselage configuration, a grid system was developed that, 
among other things, had the ability to be used for an isolated rotor configuration with minimal 
changes to the grid system. For this isolated rotor configuration, several minor changes were made 
to the grids. These changes are as follows: 

1. The nose grid has been eliminated, 

2. The tail grid has been eliminated, 

3. The collar grid has been eliminated, 

4. The top grid has been eliminated, 

5. The fuselage grid has been replaced by a volume field grid that maintains the outer boundary 
shape of the original fuselage grid. 

Since the nose, tail, collar, and top grids from items 1, 2, 3, and 4 above were “overset” into other 
background grids, it was only required to delete these grids and leave the remaining background 
grids unchanged. The only new grid in the system is the volume grid that replaces the fuselage 
grid. All other grids in the system are unchanged. In essence, what remains, is a set of volume 
grids, identical to the full configuration, minus the fuselage surface. However, now, the new rotor 
boundary condition in OVERFLOW will be applied in the rotor grid as discussed in chapter 4. 
As with the isolated fuselage configuration, the sting mount and the wind tunnel walls are not 
modeled. This new grid system is shown in figure 7.2; every second grid line has been removed 
from this figure for clarity. The farfield extent of the grid is the same as that in figure 6.3. 


7.4 Time Averaged Computation 

The isolated rotor computations will be started with a prediction of the time averaged flowfield due 
to the rotor. For this, the time averaged rotor boundary condition discussed in chapter 4 is used. 
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This is done to minimize the number of time accurate computations that must be executed. Since 
this time averaged flow field is only an intermediate stage of the computation, and since there is 
no time averaged induced inflow data for the isolated rotor test stand available, no comparisons are 
made to experimental data. 


7.5 Time Accurate Computation 

Once the time averaged computation discussed above has been obtained, OVERFLOW is restarted 
in a unsteady mode and the method is continued until a periodic solution is achieved. The state of 
“periodicity” in this context refers to the lack of differences in the induced inflow ratio predictions 
for successive blade passages. Boyd and Barnwell [24] explored the sensitivities of several param- 
eters in OVERFLOW for the time accurate computations of this type. Following the conclusions 
of that paper, the computations are executed with six Newton sub-iterations and no viscous terms. 
Since fourth order central flux differencing in space is now available in OVERFLOW, that is used 
as well. 

Figure 7.3 compares the measured and predicted unsteady induced inflow in directions parallel 
and perpendicular to the rotor tip path plane at location “J” as shown in figure 3.18. This position 
is located radially at r/R — 0.80 and azimuthally at \j/ = 84° . Figure 7.3a is the induced inflow 
comparison in the direction parallel to the tip path plane (the u-velocity, or “inplane” velocity) 
and figure 7.3b is the induced inflow comparison in the direction perpendicular to the tip path 
plane (w-velocity, or “out-of-plane” velocity). For clarity, both of the comparisons show only the 
unsteady components of induced inflow. It can be seen that the unsteady u-velocity comparison 
is excellent in magnitude, phase, and waveform shape. The w-velocity predictions match the 
phase and waveform shape well, but slightly under-predicts the magnitude of the pulse. The slight 
discrepancies in the phase for both plots may be explained by realizing that the location of the 
entire experimental signal has an error band that is equal to of the azimuthal resolution at 
which the data was acquired. This error band is equivalent to ±1.4° in azimuth. Here, the data has 
been plotted at the center of the error band. Also, the position of the measurement location in space 
and the position of the point in space used for the prediction comparisons do not exactly coincide 
since there is not a computational grid point that falls exactly on the measurement location; the 
closest point in the grid to the measurement location is shown. Comparing the coordinates for 
location J in the measurements and in the predictions, the error in radial location is less than 1% 
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of the rotor radius; the error in azimuthal location is less than 0.4° , and the vertical location error 
is negligible (-C 1% of the rotor radius). 

Even though the distributed time averaged data was not measured in this experiment, it is worth- 
while to examine this quantity from the predictions in the same manner as done by Boyd and 
Barnwell [24], Figure 7.4 shows a contour plot comparison of the time averaged induced in- 
flow (w-velocity). Figure 7.4a shows the out-of-plane component of the measured, time averaged 
induced inflow ratio for the rectangular planform rotor shown in chapter 3. This is the same exper- 
imental data from figure 3.8. As discussed in chapter 3, this experimental data contains a fuselage, 
whereas, the current prediction does not. Figure 7.4b shows the current, predicted, time aver- 
aged induced inflow ratio which as been filtered to roughly match the experimental data frequency 
content as discussed in Appendix A. 

Comparing these two figures, it can be seen that the prediction contains the same general features 
of the measured data. It can be seen that there is an upwash on the front of the disk, that there 
is a downwash at the rear of the disk with more concentrated downwash in the first and fourth 
quadrants, and that the magnitudes of the downwash are similar. These features and the level 
of prediction shown are similar to those shown in reference [24] for a different rotor system. In 
addition, figure 7.4b shows a curious feature that has been seen before in the literature for both 
the current method and for the purely GDWT [25]. This feature, seen on the retreating side of the 
rotor in the fourth quadrant near the blade root, appears to be a region where the induced inflow is 
quite small compared to the induced inflow surrounding that region; this feature is not seen in the 
measured data. That feature can be explained by examining the rotor loading in that region. Since, 
physically, the rotor induced inflow is generated as a fluid reaction to the rotor loading distribution, 
and since the rotor loading is very small in that region due to the very low dynamic pressure there, 
the induced inflow is small in that region. Chapter 8 will show that this feature is affected by the 
isolated rotor assumption made in previous investigations. 


7.6 Observations 

The inplane and out-of-plane components of unsteady induced inflow at points above, but still 
close to, the rotor plane are well predicted and are consistent with previous literature. In addi- 
tion, though time averaged results are not available for the isolated rotor configuration, the time 
averaged predicted results are consistent with previously published time averaged results at similar 
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Figure 7.1: IRTS in the NASA Langley Research Center 14- by 22-Foot Subsonic Tunnel. 
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Figure 7.3: Measured and Predicted Induced Inflow in Two Directions at Location J for the IRTS. 


David D. Boyd, Jr. 


Chapter 7. Results: Isolated Rotor 


93 


measured 



downwash 

upwash 

predicted 



Figure 7.4: Measured and Predicted Time Averaged Induced Inflow Ratio for the IRTS (Measured 
Data Includes Fuselage). 



Chapter 8 


Results: Rotor/Fuselage 


8.1 Introduction 

Chapters 6 and 7 showed various predictions made with the current computational model on the 
two separate components of a full configuration: an isolated fuselage and an isolated rotor. In 
this chapter, these two components are combined into a single unit using the entire computational 
model. For comparisons between measured and predicted quantities, several experimental data 
sets are used since there are no available data sets that contain measurements of all quantities of 
interest. 


8.2 Experimental Setup 

The experimental setup used is a combination of the IRTS discussed in chapter 7 and the 2-meter 
ROBIN fuselage [33, 34, 35, 36]. This configuration is shown in figure 8.1. In this particular 
test, no LV data were taken as was done in the IRTS test discussed previously. Instead, unsteady 
fuselage pressures were measured at several locations on the fuselage. The locations included a 
row of pressure taps near the top centerline of the fuselage. This row of taps nominally followed 
the top centerline, but were offset 0.25 inches to the advancing side of the fuselage. This offset was 
needed since the construction of the fuselage shell consisted of two halves which overlapped on 
the centerline of the fuselage. In addition, there were six unsteady pressure taps at fuselage station 
“e” (see figure 6.5) at several vertical locations on the left and right sides of the fuselage. Three of 
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these taps were on the left side and three on the right side of the fuselage. The unsteady pressure 
data was processed at 256 points per rotor revolution. 


8.3 Computational Grid System 

The grid system for the combined rotor/fuselage configuration is similar to that used in chapter 6, 
figure ??. Since the rotor grid was already included in the isolated fuselage configuration, no grid 
changes are needed to proceed directly from the isolated fuselage steady state results in chapter 6 
to the time averaged rotor/fuselage configuration. However, for the time accurate computations, 
slight modifications were made to the grid which do not affect the grid shape or distribution of 
points. These modifications were made to accommodate a solution scheme in OVERFLOW that is 
very efficient in the unsteady mode. 

For the unsteady rotor/fuselage configuration computations in OVERFLOW, the LU-SGS scheme 
is used. However, current implementation of this scheme does not allow for the use of spatially 
periodic boundary conditions. In OVERFLOW, spatially periodic grids and spatially periodic 
boundary conditions, have identical first and last planes. Figure 8.2(a) represents a schematic of 
a spatially periodic grid with k being the index in the periodic direction. It can be seen that k = 1 
and k = kmax are actually the same point; for a three dimensional grid, these would be planes. 
To account for the fact that the LU-SGS scheme cannot be used with spatially periodic boundary 
conditions, the periodic grids in the rotor/fuselage grid system (excluding the rotor grid), are con- 
verted to grids that can still pass information to themselves through an extended overlap region. 
These grids then use a new boundary condition combination, created by using several existing 
OVERFLOW boundary conditions in succession, to pass information spatially in a manner similar 
to the traditional spatially periodic boundary condition. In this new boundary condition combina- 
tion, the grids are slightly modified so that an overlap, or one-to-one match, of four grid planes 
is made in each of the computational grids’ spatially periodic direction. This overlap of four grid 
planes is used to replace the single overlap plane used in the spatially periodic boundary condition. 
A schematic of this is shown in figure 8.2(b). With these new overlapped planes, the “copy-to” 
and “copy-from” boundary conditions in OVERFLOW are used to pass information from one end 
of the computational grid to the other to simulate a spatially periodic condition. In table 8.1, the 
“copy-to” and “copy-from” columns show which information is passed between computational 
grid planes in the computational ^-direction (the spatially periodic direction used in these grids). 
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Table 8.1: LU-SGS “Periodic” Boundary Condition 


Copy-from 

Copy-to 

( k value) 

( k value) 

kmax— 3 

1 

kmax — 2 

2 

3 

kmax — 1 

4 

kmax 


There were some configuration anomalies in the experimental setup for the rotor/fuselage con- 
figuration that were also incorporated into the grid system. First, due to offsets in the wind tunnel 
floor and ceiling, the rotor center of rotation was offset approximately two inches to the advanc- 
ing side of the centerline of the fuselage. This placed the rotor rotation center at approximately 
the advancing side edge of the pylon atop the fuselage. Second, the fuselage shell was yawed 
at approximately 0.76° nose-left. Both of these anomalies have been accounted for in the grid 
system. 


8.4 Time Averaged Computation 

The isolated fuselage configuration results shown in chapter 6 are used as a starting point for the 
time averaged rotor computations with the combined configuration. These calculations involve 
using the time averaged rotor boundary condition discussed in chapter 4. As before, this portion 
of the calculation is an intermediate step between the steady state isolated fuselage computation 
and the unsteady rotor/fuselage configuration. It uses the same grid system defined above for the 
steady state computation. Since this is an intermediate stage of the computation, only the force 
convergence history will be examined to show that this portion of the computation is well behaved. 
Figure 8.3 is a continuation of figure 6.4 with the new portion beginning at iteration number 601 
and ending at iteration number 1 100. 
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8.5 Time Accurate Computation 

For the time accurate computations, OVERFLOW is executed in an unsteady mode using the LU- 
SGS scheme to invert the left hand side of the equations using the new overlapping grid scheme to 
eliminate the spatially periodic boundary conditions. The unsteady rotor model boundary condition 
is applied as discussed in chapters 4 and 7. These computations are executed until a periodic result 
is obtained in the unsteady pressures on the fuselage. 


8.6 Coupled Model Predictions 

In subsequent sections, the new, coupled, computational model predictions will be presented for 
several iterations. The first subsection will show the unsteady fuselage pressure predictions for the 
first iteration 1 vs the corresponding measured quantities from the current rotor/fuselage configura- 
tion experiment. The second subsection will compare measured and predicted, unsteady induced 
inflow ratios from the first iteration. The third subsection will compare predicted, time averaged 
induced inflow ratio contours from the first iteration, which have been derived from the unsteady 
model, to the corresponding measured induced inflow ratio contours, taken from the experimental 
data presented in chapter 3. The fourth subsection will compare measured and predicted lateral and 
longitudinal induced inflow ratios. Subsequent sections and subsections will show the comparisons 
for the second and third iterations, including the coupling. Finally, an examination of the pressure 
contours on the fuselage surface will be made, as well as an examination of the convergence of the 
method in terms of the trim pitch settings as a function of iteration. 


8.6.1 Iteration 1 

Unsteady Fuselage Pressures 

Figure 8.4 shows a comparison between the measured and predicted unsteady modified pressure 
coefficient on the top centerline of the ROBIN fuselage for various downstream locations (x- 
direction). Plotted on the vertical axis is the negative of the modified pressure coefficient, d p . 

^ere, an iteration refers to one pass through the rotor loading model followed by one pass through the ro- 
tor/fuselage flowfield model. Typically, only two rotor revolutions are required in each rotor/fuselage flowfield model 
iteration. 
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Along the horizontal axis is the reference blade location in the rotor azimuthal coordinates. Since 
the measurement location is fixed in the non-rotating frame, this horizontal axis can also be viewed 
as a temporal axis spanning one rotor revolution. The modified pressure coefficient used in these 
comparisons is defined by the following: 


, ioo(p-p..) 
p " ip m 2 


( 8 . 1 ) 


where the typical velocity in the denominator of the Cp definition had been changed to OR. This 
modified C p is designated d p . Using a superscript u simply denotes that this is the unsteady compo- 
nent. This change in definition is needed since, for a rotorcraft, a freestream velocity approaching 
zero is possible as the hover condition is approached. If the normal definition of C p were used, 
the values of C p would approach infinity as the hover condition is approached. The factor of 100 
on the d p simply serves as a convenient factor to scale the modified pressure coefficient function. 
Comparing the standard C p definition and the d p definition above, the following relation holds 
between C p and d p : 


C p =l00/jlC p (8.2) 

where p„ is the freestream advance ratio. To show the relation between these quantities more 
clearly, consider a perfect fluid. For a perfect fluid, stagnation occurs at — C p = —1.0. Using the 
current flight condition for the rotor/fuselage configuration, and using equation (8.2), this equates 
to -d p = -5.29. Referring this to figure 8.4, a value of -d p = -1.0 corresponds to - C p « -0. 19. 

Examining features in figure 8.4, it can be seen in that there is a dominant blade passage event at 
a frequency of four pulses per rotor revolution in the measured pressure signature. Here, the phase 
of the pressure signal is well matched, and the magnitude is slightly over-predicted. 

Figure 8.5 shows the measured and predicted -d p as a function of the reference blade location. 
Whereas figure 8.4 showed the predictions on the top centerline of the fuselage, figure 8.5 shows 
the predictions on the retreating and advancing sides of the fuselage ( i.e ., the left and right sides) for 
various vertical locations at the same downstream location. It can be seen that the magnitude and 
phase are predicted well for the advancing side locations; however, the retreating side amplitudes 
are slightly over-predicted. 
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Time Accurate Induced Inflow Ratios 

Figure 8.6 shows the in-plane and out-of-plane components of unsteady induced inflow for the 
same location as that in figure 7.3. Figure 8.6, however, is a comparison between the unsteady 
induced inflow velocities for an isolated rotor (obtained from figure 7.3) and that for the first 
iteration of the method, including the fuselage. It can be seen that the presence of the fuselage has 
little effect on the unsteady induced inflow components at these particular locations. 


Time Averaged Induced Inflow Ratios 

At this point, it is worthwhile to check the validity of the time averaged results, obtained by time 
averaging the time accurate calculations, by comparing them to measured data. 

Figure 8.7 shows a comparison between the measured and predicted induced inflow ratio per- 
pendicular to the rotor tip path plane. For each of the predictions, the same pressure distribution is 
used in the rotor boundary condition. Figures 8.7b and 8.7c show the predicted induced inflow ra- 
tio for the isolated rotor and the predicted induced inflow ratio for the rotor/fuselage combination, 
using the same unsteady rotor boundary condition and pressure information. Several improve- 
ments can be seen in the induced inflow ratio prediction between the isolated rotor prediction and 
the rotor/fuselage combination. First, the inflow anomaly, seen in the isolated rotor prediction in 
the fourth quadrant near the blade root, is not present in the rotor/fuselage combination predic- 
tion. The induced inflow in this region now matches the measured data in that region. Second, 
the induced inflow near the forward portion of the rotor disk matches the experimental data well. 
For example, examining the contour line of “zero” induced inflow in the measured data and in the 
rotor/fuselage combination prediction, it can be seen that the upwash on the forward section of the 
disk now extends toward the center of the rotor. This upwash region is generated by the presence 
of the fuselage and is therefore not predicted in the isolated rotor prediction. 

Figure 8.7d shows the difference between the induced inflow ratio for the rotor/fuselage combi- 
nation and that of the isolated rotor. This difference shows the effects of the fuselage on the time 
averaged flowfield. As expected, there is an increased upwash near the nose of the fuselage and 
forward section of the pylon and an increased downwash behind the pylon. The regions of addi- 
tional induced inflow in this figure show the source of the improved induced inflow predictions 
shown in figure 8.7c. 

Figure 8.8 shows a comparison between the measured and predicted induced inflow ratio parallel 
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to the rotor tip path plane for the same conditions shown above. Figures 8.8b and 8.8c show the 
predicted induced inflow ratio parallel to the tip path plane for the isolated rotor and the predicted 
induced inflow ratio parallel to the tip path plane for the rotor/fuselage combination. Comparing 
the figures, it can be seen that the rotor/fuselage combination better predicts the induced inflow 
ratio. For example, examining the induced inflow contour line with a value of 0.015, it can be 
seen that the isolated rotor prediction is almost symmetric left-to-right, whereas the measured and 
predicted induced inflow ratios are both skewed at about 45° counterclockwise to the oncoming 
flow. Here, again, it can be seen that inclusion of the fuselage is necessary to correctly predict the 
time averaged induced inflow ratio. 

Figure 8.8d shows the difference between the induced inflow ratio parallel to the tip path plane 
for the rotor/fuselage combination and that of the isolated rotor. This difference, again, shows the 
effect of the fuselage on the time averaged flowfield and the sources of the improved predictions 
with the fuselage included. As expected, there is a decrease in u-velocity near the nose of the 
fuselage and an acceleration near the rear of the fuselage pylon. It can be seen that, to capture the 
u-velocity distribution accurately, the fuselage must be included in the computation. 


Lateral/Longitudinal Induced Inflow Ratios 

Figure 8.9 compares the lateral and longitudinal subsets of the measured and predicted induced 
inflow ratio. The measured data comes from figure 8.7a, while the predicted data is taken from 
the combined rotor/fuselage case in 8.7c. Comparing the prediction here to the prediction using 8 
harmonics in figure 3.9, it can be seen that the longitudinal induced Inflow ratio is well predicted 
using the combined rotor/fuselage model. Whereas the peak-to-peak amplitude of the lateral in- 
duced inflow ratio was over-predicted by the GDWT isolated rotor model, it is well predicted for 
the combined rotor/fuselage model. However, there is a slight over-prediction of the inflow ratio 
near the advancing and retreating blade tip regions. 


8.6.2 Iteration 2 

Unsteady Fuselage Pressures 

With the first iteration complete, induced inflow corrections are computed as described in chapter 
5. With these induced inflow corrections, the GDWT is used to recompute the rotor loading. 
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With this new rotor loading, OVERFLOW is used to recompute the rotor/fuselage flowfield. This 
computation, or “iteration 2”, is presented in the current and subsequent sections. 

Figure 8.10 shows the top centerline modified pressure coefficient in a manner similar to figure 
8.4. Comparison of these two figures shows only very small differences. 

Figure 8.1 1 shows the the modified pressure coefficient on the retreating and advancing sides of 
the fuselage in a manner similar to figure 8.5. Comparison of these two figures, again, shows only 
very small differences. 

Time Accurate Induced Inflow Ratios 

Figure 8.12 shows the in-plane and out-of-plane components of unsteady induced inflow for the 
same location as that in figure 7.3. As was seen in the first iteration, the fuselage has little influence 
on the unsteady components of induced inflow at this point in the flowfield. 


Time Averaged Induced Inflow Ratios 

Figure 8.13 shows the same information as figure 8.7, but this is the second iteration. Figures 8.13a 
and 8.13b show the identical information as that in figures 8.7a and 8.7b. Plots 8.13c and 8.13d 
are the new parts of this figure. Comparing figures 8. 13c and 8.7c, it can be observed that there are 
only small changes between these figures except on the advancing side of the rotor disk in the first 
quadrant. In the first iteration, the time averaged induced inflow in that region was over-predicted. 
In the second iteration, it is seen that the time averaged induced inflow quantities are now closer 
to the measured values. However, overall, the shape of the contours did not change significantly. 
Comparing figures 8.13d and 8.7d, it can be observed that there are only small changes between 
these difference plots. Thus the effect of the fuselage on the time averaged induced inflow is quite 
similar between the first and second iterations. 

Figure 8.14 shows the same information as figure 8.8, except it is for the second iteration. As 
before, figures 8.14a and 8.14b show the same information as figures 8.8a and 8.8b. Comparing 
figures 8.14c and 8.8c, it can be seen, again that there are only small differences between the first 
and second iterations. This time, however, the small differences that do occur are on the advancing 
side in the first part of the second quadrant of the rotor disk. A similar conclusion holds for figures 
figures 8.14d and 8.8d. 
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Lateral/Longitudinal Induced Inflow Ratios 

Figure 8.15 compares the lateral and longitudinal subsets of the measured and predicted induced 
inflow ratio in a manner similar to figure 8.9. The measured data comes from figure 8.13a (and is 
the same as in figure 8.7a), while the predicted data is taken from the combined rotor/fuselage case 
in 8.13c. The predictions shown in this figure are similar to those shown in figure 8.9, and similar 
trends are present. 

Even though differences are small between the first and second iterations, a third iteration is also 
presented to show that the iteration process has converged. 

8.6.3 Iteration 3 

Unsteady Fuselage Pressures 

After iterating once again through the coupling technique, the rotor loading model, and the ro- 
tor/fuselage flowfield model, the third iteration is complete. Figures 8.16 and 8.17 once again 
show the unsteady modified pressure coefficient on both the top centerline and the sides of the 
fuselage. As with the comparison between the first and second iterations, there is no significant 
difference between this iteration and the previous iteration. 


Time Accurate Induced Inflow Ratios 

Again comparing two components of unsteady induced inflow ratio for the current iteration to the 
isolated rotor results in figure 8.18 shows no significant difference between the two results. In 
addition, there are no significant differences between this iteration and the previous iteration. 


Time Averaged Induced Inflow Ratios 

Examining the in-plane and out-of-plane time averaged induced inflow ratios in figures 8.19 and 
8.20, as was done for the second iteration, it can be seen that there are no significant differences 
between the current iteration and the previous iteration. 

Since there are no significant differences between the second and third iteration in the unsteady 
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modified pressure coefficient, in the time accurate inflow ratio, or in time averaged induced inflow 
ratio, it can be concluded that the iteration process is complete. 

Lateral/Longitudinal Induced Inflow Ratios 

Figure 8.21 compares the lateral and longitudinal subsets of the measured and predicted induced 
inflow ratio in a manner similar to figures 8.9 and 8.15. The measured data comes from figure 
8.19a (and is the same as in figures 8.7a and 8.13a), while the predicted data is taken from the 
combined rotor/fuselage case in 8.19c. Again, the predictions shown in this figure are similar to 
those shown in figures 8.9 and 8.15, and similar trends are present. 

8.7 Examination of Pressure Contours 

In the preceding sections, a detailed examination of the unsteady component of the modified pres- 
sure coefficient, of the time averaged inflow ratio, and of the unsteady inflow ratio was presented 
for several iterations of the current method. For the final result (at the end of the third iteration), 
it is constructive to examine pressure distribution on the fuselage with and without the effects of 
the rotor. Figure 8.22(a) is a contour plot of surface pressure coefficient on the isolated fuselage 
configuration presented in chapter 6. It can be seen there is a typical stagnation region on the nose 
of the fuselage. 

Figure 8.22(b) is a contour plot of the time averaged surface pressure coefficient for the ro- 
tor/fuselage configuration. This plots contains data from the “iteration 3” above. As with figure 
8.22(a), there is a stagnation region on the fuselage nose region and behind the pylon. However, 
comparing figures (a) and (b), it can be seen that the gross time averaged effect of the rotor is to 
increase the pressure coefficient on the sides of the fuselage and downstream of the pylon. Figures 
8.23 and 8.24 are top views of figures 8.22(a) and (b), respectively. 

Figure 8.25 shows a set of surface pressure coefficient contour plots similar to those presented 
in figure 8.22. For comparison purposes, figure 8.25(a) shows the same data as presented in figure 
8.22(b). Figure 8.25(b) shows the surface pressure coefficient for the instant at which the reference 
blade is at the \\f = 0° location. All four blades are plotted as rectangles extending from the blade 
root at r/R — 0.24 to the blade tip at r/R= 1.0 and are drawn “to-scale” in their actual position and 
orientation with respect to the fuselage; The blade chord is drawn “to-scale” as well. In this contour 
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plot from the unsteady computation, a large pressure pulse can be seen on the upper surface of the 
tail boom. This pressure pulse is a direct result of the rotor blade at the \|i = 0° location passing 
over the tail boom. Figure 8.26 shows this pressure pulse in a top view of figure 8.25(b). 

Figures 8.27 to 8.35 present surface pressure coefficient contour plots at azimuth locations of 
\j r = 15° ,30° ,45° ,60° ,75° , and 90° . These plots again show a pressure pulse traveling on and 
around the fuselage. It can be seen that this pulse motion is correlated with the motion of the 
blade. This point is made especially clear in the “Top View” plots presented here. Examination of 
this set of plots reveals that, locally, the unsteady surface pressure coefficient can be substantially 
higher than the time averaged surface pressure coefficient and that the pulses are typically of short 
temporal duration. In addition, the asymmetry of the unsteady loading can be seen clearly in the 
“Top View” figures. 


8.8 Iteration Effects on Rotor Trim 

The effect of the coupled computation on the rotor trim can be assessed by examining the pitch 
control settings required at each stage of the iteration at the end of the trim procedure in the rotor 
loading computation. These blade pitch settings are a function of rotor azimuth, are referenced to 
the 0.75R radial location on the blade, and are defined as follows: 

0(\|/) = 0o + 9 c cosv)r+0iSin\)/ (8.3) 

where 0o is the collective pitch, 0 C is the longitudinal pitch, and 0 5 is the lateral pitch. The first 
row in table 8.2 shows the measured collective, longitudinal, and lateral pitch settings, along with 
the magnitude and phase of the longitudinal/lateral combination. All quantities in table 8.2 are in 
degrees. The magnitude of the first harmonic of pitch, 0j, and phase of the first harmonic of pitch, 
\jq can be defined as follows: 


0(v) 

= 0o + 0lCOs(\}/-\l/i) 

(8.4) 

01 

= \/ e ? +6 ‘ 2 

/ f\ \ 

(8.5) 

Yi 

= 

(8.6) 
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where \|ti is placed in the range 0° < vjq < 360° . 

Table 8.2: Pitch Control Settings [deg] as a Function of Iteration Number 


Type 

Iteration 

e 0 

0c 

0, 

01 

Vi 

Measured 

N/A 

8.16 

1.52 

-4.13 

4.40 

290.2 

Isolated Rotor 

1 

6.84 

1.45 

-3.42 

3.71 

293.0 

Rotor+Fuselage 

2 

8.06 

2.73 

-3.69 

4.59 

306.5 

Rotor+Fuselage 

3 

8.02 

2.67 

-3.71 

4.57 

305.7 


From this table, it can be seen that the collective pitch setting matches the measured collective 
pitch for the rotor/fuselage combination cases, but is under-predicted for the isolated rotor case as 
expected. It can also be seen that the magnitude of the first harmonic of pitch is well matched for 
all of the iterations, especially for the iterations that include the fuselage in the computation. A 
plot of these quantities vs iteration number are shown in figure 8.36. From this figure, it can be 
seen that these settings have converged in just two iterations. For the computations that include 
the complete configuration, a phase difference of approximately 16° can be seen. The explanation 
of this phase difference can be attributed to the assumption that the blade flap hinge offset is at 
the center of rotation. To show this, it is only necessary to examine the effect of the blade hinge 
offset on the rotor flap response. Gessow and Myers [54] showed that the rigid flapping natural 
frequency of a rotor is a function of blade hinge offset and is given by the following formula: 


co„ = 



(8.7) 


where co„ is the flapping natural frequency in cycles per revolution and h is the hinge offset location 
as a fraction of rotor radius. Examination of equation (8.7) for h — 0 shows that the flap natural 
frequency co„ = 1.0 per revolution and for h = 0.06, the flap natural frequency is co„ = 1.044 per 
revolution. The difference in these two natural frequencies is A(0„ = 0.044 per revolution. Since 
there are 360° in one revolution, the natural frequency difference equates to approximately 16° 
of rotor azimuth. This means that in the experiment, which included the blade hinge offset, the 
phasing of the first harmonic of pitch should preceed the predicted phasing by approximately 16° 
of rotor azimuth. This is precisely the amount seen in table 8.2. 
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8.9 Resource Usage Summary 

Table 8.3 lists the resource usage for each component of the current rotor/fuselage combination 
computations using the GDWT with 8 harmonics and 128 azimuth steps per revolution. The items 
listed are the CPU times in hours and minutes [hr: min] required for each iteration stage, main 
memory required in mega-words [mw], and which machine was used for each phase. 


Table 8.3: Resource Usage as a Function of Iteration Number 


Type 

CPU [hr:min] 

Memory [mw] 

Machine 

GDWT (with 8 harmonics, 128 azimuths) 

<0:07 

5 

Cray C-90 

Isolated Fuselage 

3:43 

18 

Cray C-90 

Time Averaged, Isolated Rotor 

2:16 

15 

Cray C-90 

Time Accurate, Isolated Rotor (per rev.) 

1:52 

17 

Cray C-90 

Time Averaged, Rotor/Fuselage 

3:38 

18 

Cray C-90 

Time Accurate, Rotor/Fuselage (per rev.) 

2:42 

22 

Cray C-90 


Each of the rows in table 8.3 is for a particular component of the method for this particular case. 
Since execution times will differ for each particular case, these should only be used as reference 
quantities. It should be noted that OVERFLOW does not have a convergence criterion that halts 
execution at a particular convergence level. It is executed for a specified number of iterations or 
time steps. Each of the components involving time accurate computations are given as the CPU 
time required to execute each rotor revolution. For the case presented here, each time accurate 
stage of the computation was executed for two complete rotor revolutions to assure periodicity, 
even though periodicity may be achieved in fewer actual blade passage events. Thus, the CPU 
times above are conservative values for this case. Also, note that the GDWT is executed in a 
conservative manner at the same resolution required by OVERFLOW. Thus, CPU time for the 
GDWT computations is conservative as well. To compute the total execution time for the isolated 
rotor results and the first iteration of the rotor/fuselage configuration, the following formulae can 


be used: 



w 








w 
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Isolated Rotor: 

CPU Time 

+ 

+ 


GDWT 

(Time Averaged, Isolated Rotor) 
2(Time Accurate, Isolated Rotor) 
6:07 (Cray C-90) 


Rotor/Fuselage, Iteration 1 : 
CPU Time 

+ 

+ 

+ 


(Isolated Fuselage) 

GDWT 

(Time Averaged, Rotor/Fuselage) 
2(Time Accurate, Rotor/Fuselage) 
12:52 (Cray C-90) 


where the factor of two on the time accurate components indicates that two complete revolutions 
were executed for this particular case. 


8.10 Observations 

This chapter has presented results from the full hybrid method for a rotor/fuselage configuration. 

From the evidence presented above, several conclusions can be drawn as follow: 

• The unsteady modified pressure coefficient on the top centerline of the fuselage and on the 
sides of the fuselage shown above, are insensitive to the small trim changes made by the time 
averaged induced inflow corrections which account for the presence of the fuselage. 

• The in-plane and out-of-plane unsteady induced inflow velocity components are also insen- 
sitive to the small trim changes made by the time averaged induced inflow corrections which 
account for the presence of the fuselage. 

• The time averaged induced inflow velocities are improved through the iteration/coupling 
process presented above. For the case presented here, it appears that one iteration is a good 
approximation to the correct solution and that two iterations is sufficient to capture the time 
averaged and time accurate induced inflow effects. 

• The primary effect of the fuselage on the trimmed blade pitch settings is on the collective 


v 
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pitch (for this particular case). The presence of the fuselage improves the pitch setting pre- 
dictions over those from the isolated rotor case. 

• For this flight condition, the primary effects of the rotor on the fuselage are a higher time 
averaged surface pressure coefficient below and downstream of the rotor and short duration 
surface pressure pulses imposed by the individual blade passages over the fuselage surfaces. 
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Figure 8.1: IRTS/Fuselage Configuration in the NASA Langley Research Center 14- by 22-Foot 
Subsonic Tunnel. 





V 






David D. Boyd, Jr. 


Chapter 8. Results: Rotor/Fuselage 


110 


kmax-3 "V 

kmax-2*^--*_ 

kmax-1 


\ \ 

kmax-6 
km ax-5 

kmax-4 


kmax-3 


3 y kmax 

r ^ kmax-1 
kmax-2 ( b ) 


Figure 8.2: Schematic of a Periodic Grid and Replacement for Periodic Grid in LU-SGS Scheme. 
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Figure 8.3: Lift, Drag, and Sideward Direction Force Coefficients, Time Averaged Computation. 
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Figure 8.4: Measured and Predicted Unsteady Modified Pressure Coefficient on the Top Centerline 
of the ROBIN Fuselage, Iteration 1. 
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Figure 8.5: Measured and Predicted Unsteady Modified Pressure Coefficient on the Retreating and 
Advancing Sides of the ROBIN Fuselage, Iteration 1. 
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Figure 8.6: Measured and Predicted Induced Inflow in Two Directions for an Isolate Rotor and a 
Rotor/Fuselage Combination, Iteration 1. 
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Figure 8.7: Measured and Predicted Time Averaged Induced Inflow Ratio from Time Accurate 
Computations, Iteration L 
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Figure 8.8: Measured and Predicted Time Averaged Parallel Induced Inflow Ratio from Time 
Accurate Computations, Iteration 1. 
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Figure 8.9: Measured and Predicted Lateral and Longitudinal Time Averaged Induced Inflow Ratio 
from Time Accurate Computations, Iteration 1 
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Figure 8.12: Measured and Predicted Induced Inflow in Two Directions for an Isolate Rotor and a 
Rotor/Fuselage Combination, Iteration 2. 
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Figure 8.13: Measured and Predicted Time Averaged Induced Inflow Ratio from Time Accurate 
Computations, Iteration 2. 
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Figure 8.14: Measured and Predicted Time Averaged Parallel Induced Inflow Ratio from Time 
Accurate Computations, Iteration 2. 
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Figure 8.15: Measured and Predicted Lateral and Longitudinal Time Averaged Induced Inflow 
Ratio from Time Accurate Computations, Iteration 2 
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Figure 8.16: Measured and Predicted Unsteady Modified Pressure Coefficient on the Top Center- 
line of the ROBIN Fuselage, Iteration 3. 



David D. Boyd, Jr 


Chapter 8. Results: Rotor/Fuselage 


125 


measured 

predicted 


0 90 180 270 360 

Ref. Blade Loc. (°) 


3 

2 

■ C P 0 

-1 

-2 

Retreating Side 

z=-0.07 3 

: 2 

: 1 

Retreating Side 

z=0.06 3 ' 

: 2 : 

: 1 : 

-1 

: -2 

: -1 : 

-2 : 

3 ( 

90 180 270 360 w ( 

90 180 270 360 0 

3 

Advancing Side 3 

. Advancing Side 3 


GO 

O 

■ 

0 

1 

il 

N 

z=0.06 : 

2 

2 

P 2 • 

1 

1 

: 1 ■ 

C p 0 

0 


-1 

-1 

-i * 

-2 

-2 

-2 - 

-3 

, i , . i . . ... i , , , .-i _q 



Retreating Side 
z=0.12 




90 180 270 360 

Advancing Side 
z=0.12 




0 90 180 270 360 

Ref. Blade Loc. (°) 


0 90 180 270 360 

Ref. Blade Loc. (°) 


Figure 8.17: Measured and Predicted Unsteady Modified Pressure Coefficient on the Retreating 
and Advancing Sides of the ROBIN Fuselage, Iteration 3. 
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Figure 8.18: Measured and Predicted Induced Inflow in Two Directions for an Isolate Rotor and a 
Rotor/Fuselage Combination, Iteration 3. 
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Figure 8.19: Measured and Predicted Time Averaged Induced Inflow Ratio from Time Accurate 
Computations, Iteration 3. 
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Figure 8.20: Measured and Predicted Time Averaged Parallel Induced Inflow Ratio from Time 
Accurate Computations, Iteration 3. 
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Figure 8.21: Measured and Predicted Lateral and Longitudinal Time Averaged Induced Inflow 
Ratio from Time Accurate Computations, Iteration 3 



David D. Boyd, Jr. 


Chapter 8. Results: Rotor/Fuselage 


130 



Figure 8.22: Isolated Fuselage and Rotor/Fuselage, Time Averaged Surface Pressure Coefficients 
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Figure 8.23: Surface Pressure Coefficient on Isolated Fuselage Configuration, Top View 
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Figure 8.24: Time Averaged Surface Pressure Coefficient on Complete Rotor/Fuselage Configura- 
tion, Top View 
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Figure 8.26: Time Accurate Surface Pressure Coefficient on Complete Rotor/Fuselage Configura- 
tion, \|/ = 0° , Top View 
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Figure 8.27: Time Accurate Surface Pressure Coefficients at y = 15° and 30° 
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Figure 8.28: Time Accurate Surface Pressure Coefficient on Complete Rotor/Fuselage Configura- 
tion, \jr = 15° , Top View 
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Figure 8.29: Time Accurate Surface Pressure Coefficient on Complete Rotor/Fuselage Configura- 
tion, vjr = 30° , Top View 
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Figure 8.32: Time Accurate Surface Pressure Coefficient on Complete Rotor/Fuselage Configura- 
tion, \g = 60° , Top View 
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Figure 8.33: Time Accurate Surface Pressure Coefficients at \|/ = 75° and \j/ = 90° 
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Figure 8.35: Time Accurate Surface Pressure Coefficient on Complete Rotor/Fuselage Configura- 
tion, \|/ = 90° , Top View 
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Figure 8.36: Pitch Settings as a Function of Iteration Number 



Chapter 9 


Summary 


In this research, an efficient and accurate, hybrid method, combining a rotor loading model, a 
rotor/fuselage fiowfield model, and a coupling technique, has been developed. This method uses 
the GDWT for the rotor loading model, OVERFLOW for the rotor/fuselage fiowfield model, and 
a new coupling technique developed in this research. The motivations for development of such a 
model have been discussed and the relationships between the current model and other models in 
use today have been discussed. 

A description of the Rotor Loading Model used for the first component of the model has been 
presented. In addition, the solution procedures employed in this model have been discussed, along 
with a validation study of the Rotor Loading Model. In this study, it was shown that predicted 
time averaged inflow distributions over a rotor disk matched measured time averaged inflow distri- 
butions well for a several flight conditions and for several rotor configurations. The sensitivity of 
the results to parameters (such as the number of azimuth steps and the number of harmonics used) 
was shown. In addition, predicted time accurate inflow quantities were compared to measured 
quantities. These comparisons showed that the predictions of the unsteady quantities do not match 
the experimental values well in waveform, but are comparable in magnitude and phase. It is also 
noted here that these predictions match well with previously published predictions in the literature 
for these rotor configurations. 

A discussion of the method used for the second component of the current model, the Ro- 
tor/Fuselage Fiowfield Model was also provided. The method used to solve the Navier-Stokes equa- 
tions was discussed; theory and implementation of the new time averaged and time accurate bound- 
ary conditions is discussed as well. These new rotor boundary conditions are a unique feature of the 
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current model, and this is the first time this type of condition has been used for such a computation. 

A coupling model was developed, which links the Rotor Loading Model and the Rotor/Fuselage 
Flowfield Models through the time averaged induced inflow of the rotor. This coupling model is a 
unique feature of the current model. The derivation of a novel inflow filtering technique has been 
presented also. 

Computations which use the model presented in the first several chapters are compared with 
experimental data. These comparisons include results for an isolated fuselage, an isolated rotor, 
and a rotor/fuselage configuration. Predictions of the pressure coefficient on the surface of the 
isolated fuselage were shown to match experimental data well. Predictions of the time average and 
time accurate inflow above the rotor tip path plane for the isolated rotor configuration were shown 
also to match experimental data well. 

Predictions of time averaged and time accurate inflow velocities above the rotor disk and predic- 
tions of unsteady pressure coefficients on the top centerline of the fuselage were shown to match 
experimental data well. A unique feature of this model is that predictions of the pressure coefficient 
on the sides of the fuselage match well with measured data. This is a significant accomplishment; 
previous methods have been unable to match unsteady pressures on the sides of the fuselage. 
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Appendix A 


Filtering Operation 


A consistent filtering operation can be derived from the GDWT. This appendix gives a derivation 
of the filtering operation, derived using various functions introduced in the GDWT development 
of reference [28], This operation takes a quantity, A(r,y), which is a function of the radial and 
azimuthal coordinate, and derives the expressions for coefficients of an infinite series given below. 
The A(r,y) quantity is then filtered by truncating the infinite series to a finite number of harmonics 
and shape functions. 

Starting with a given quantity that is a function of the radial and azimuthal coordinates, one 
can express that quantity as a double summation over the harmonics and shape functions in the 
GDWT as follows: 


A(r,y) = XS ft? (?) K cos (my) + &™sin(my)] (A.l) 

m n 

where A is the given function, m is the harmonic number (m = 0, 1,2. . .), n is the shape function 
number (n = m+ l,m + 3, . . .)> r and y are the radial and azimuthal coordinates, <j)™ is a function 
of r from the GDWT , and a ™ and b% are the unknown coefficients. These coefficients could be, 
in general, a function of time. In that case, the process outlined below would be followed at each 
discrete time step of interest. 

Multiply equation (A.l) by cos(py) and integrate over y from 0 to 2tt to get the following: 
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r 271 

/ A(r,\|/)cos(pv)cty 
Jo 


271 


SSC( F )[ a ^ I cos(my)cos(p\\f)dy 

m n 0 

271 

+6™ / sin(m\|/) cos(p\|/) d\\f } 
o 


(A.2) 


Using the orthogonality relations for sine and cosine functions, the right hand side of equation 
(A.2) can be rewritten. The resulting form of equation (A. 2) can be written as follows: 


00 00 

/ A(r, \|/) cos(/n|/) d\\f = 

JO m n 

= £m&(p) < A - 3 > 

n 

where m = 0, 1 , 2, for both equations above, and n = m + 1 , m + 3 , . . . , °°, for the first equation 
of (A. 3) and n ^ p + l,p + 3,. . . for the second equation of (A. 3). Also, 5 m p is the Dirac delta 
function. The function c(p) is given as follows: 


c(p) = 


n for p± 0 
2n for p = 0 


(A.4) 


Multiplying equation (A.3) by [<j)£(r) ■ r ■ a/1 — r 2 ] and integrating from r — 0 to r — 1 gives the 
following: 


Integral A 


I'^rv/l-f 2 A(r, \\f) cos(p\j/) d\y = £ J Q \ - r 1 dr 


a p n c{p) 

(A. 5) 


Into “Integral A”, substitute the definition of <j> in terms of Legendre functions, and perform a 
change of variables from r to v using the fact that v = \J 1 — r 2 . The resulting Integral A compo- 
nent of the above equations becomes: 
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J Q &(r)$P{r)fVl-r 2 dr 

- fT [J-fo 

~ l v v 4 V«fV"? 




4 s/h£h% 


J mp 


where H ™ is defined in the GDWT as follows: 


K = 


(n + m— 1 )!!(n — m— 1)!! 


(n + m ) ! ! (n — m ) ! ! 
where the double factorial is defined in Peters, et al. [26], as follows: 


(n)(n - 2)(n — 4) . . 

•(2) 

/<w 

n —even 

(n)(n-2)(n-4).. 

•( 1 ) 

/or 

n —odd 

1 


f or 

n —0 

1 


for 

« =— 1 

, -1 


for 

n =— 3 


Substituting equation (A.6) back into equation (A.5) gives the following: 


KJ 


[ <jj £(r)n/ 1 - r 2 [ A(r,\\f)cos(p\\f)d\\f dr 

Jo H Jo 


= I 

n 


nc(p) p 

4HS a Q 


a p n c{p) 


Solving equation (A. 8) for the unknown coefficient a p one gets: 


2 P = _ll±_ 

q nc{p) 


4 H p r 


j §q{r)r\J 1 — r 2 jf A(r,\|/)cos(p\|/) d\\l 


(A.6) 


(A. 7) 


(A.8) 


dr 


(A.9) 
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noting that p = 0, 1,2, . and q = p + l,p + 3, .. Equation (A.9) gives the unknown a p q 
coefficients of equation (A.l) given a function A(r,\|/). A similar derivation can be performed for 
the bq coefficients, which results in the following: 


b r = 12L 
* n Up) 


r\ _ f r2n 

/ €(rWl-r 2 / A(r,\y)sin(/?\|/)d\|/ 

Jo H L-/0 


dr 


(A. 10) 


noting that p = 1 , 2, . . . , °o and q = p + 1 , p + 3 , . . . , For this application, the integrals in equations 
(A.9) and (A. 10) are computed using the trapezoidal rule. 

In order to filter the function to a contain a particular frequency content, it is necessary to 
truncate the summation to limit the number of a q and b q coefficients to a specified number of 
harmonics such that p = 0,1,2 ,...,Pmax for the a q coefficient and p = 1,2 ,...,p max for the b q 
coefficient. In this implementation, the above limits set the number of shape functions to q = 
p+ l,p + 3,...,Pmax ■ Renaming the indices, and applying truncation to the summations, the 
filtered value of A(r,\(/), shown below as A(r,\|/), is as follows: 


m ma x m max^~ l 

A(r,\|/)=£ £ C(f) Kcos(my)+ ^sin(mv)] (A.ll) 

m n 

where, m is the harmonic number index, n is the shape function index, and m max is the number of 

harmonics used. It should be noted that the number of shape functions used is equal the number of 

harmonics used, plus one ( n max = m max + 1)- Thus, compacting the above notation and referring to 
the terminology of chapter 5, the following filtering operation is defined: 
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Abstract 

A novel unsteady rotor-fuselage interactional aerody- 
namics model has been developed. This model loosely 
couples a Generalized Dynamic Wake Theory (GDWT) 
to a thin-layer Navier-Stokes solution procedure. This 
coupling is achieved using an unsteady pressure jump 
boundary condition in the Navier-Stokes model. The 
new unsteady pressure jump boundary condition mod- 
els each rotor blade as a moving pressure jump which 
travels around the rotor azimuth and is applied between 
two adjacent planes in a cylindrical, non-rotating grid. 
Comparisons are made between measured and predicted 
time-averaged and time-accurate rotor inflow ratios. Ad- 
ditional comparisons are made between measured and 
predicted unsteady surface pressures on the top center- 
line and sides of the fuselage. 

Introduction 

It is well known that rotorcraft aerodynamics is a com- 
plicated topic. Due to the combination of various sys- 
tems associated with rotorcraft, these aerodynamic phe- 
nomena are unsteady, even in level, unaccelerated flight. 
Complicating these issues are the facts that typical rotor- 
craft in service today have bluff aft regions, which can 
lead to large regions of flow separation, and that there can 
be significant aerodynamic interaction or interference be- 
tween the rotating and non-rotating components of the 
system. 
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When designing a new rotorcraft, as with any flight 
vehicle, an understanding of the aerodynamic environ- 
ment, including aerodynamic interaction of the different 
vehicle components, is essential. These interactional ef- 
fects have been known and categorized for many years. 
This paper focuses on the “rotor-fuselage” and “fuselage- 
rotor” subsets of the categories offered by Sheridan and 
Smith. 1 In practice, information on specific interactional 
effects may be obtained using any combination of wind 
tunnel testing and/or computational modeling. 

Wind tunnel testing has been relied upon heavily in 
designing new rotorcraft and diagnosing and correcting 
aerodynamic anomalies discovered on actual flight vehi- 
cles because computational modeling of rotorcraft aero- 
dynamics is still in its infancy and lags well behind the 
computational capabilities used for fixed wing vehicle 
modeling. Several factors have led to this situation. One 
of these is the fact that, as mentioned above, even in 
level, unaccelerated flight, a rotorcraft is operating in an 
unsteady aerodynamic environment due to the rotation 
of the rotor system. A fixed wing aircraft in the same 
situation would be in a steady state environment. The 
computational implication of this is that a complete ro- 
torcraft simulation would necessarily be a time-accurate 
computation, whereas the fixed wing simulation could 
be a steady-state computation. Another factor is asso- 
ciated with the vastly different time and length scales 
associated with rotorcraft. Some unsteady aerodynamic 
events, such as blade-vortex interaction, occur at length 
scales that are a small fraction of a blade chord and at 
time scales that are equivalent to a tiny fraction of a rotor 
revolution. To capture these effects, very small time steps 
would be required. However, determining the trim state 
of a rotorcraft requires balancing the gross forces on the 
rotorcraft that have a length scale on the order of the ro- 
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Figure 1 . Analysis types for coupled solutions. 

tor radius (i.e., many chord lengths) over a relatively long 
time scale, equivalent to a number of rotor revolutions. 
The computational implication of these vastly different 
time scales is that a time-accurate simulation would need 
to be executed for many time steps. 

There are a number of methods available for compu- 
tation of the interactional aerodynamic effects associated 
with rotorcraft. Figure 1 categorizes these methods into 
three areas: Singularity Methods, Hybrid Methods, and 
Computational Fluid Dynamic (CFD) Methods. Each of 
these methods has been used in the past for computa- 
tion of rotorcraft interactional aerodynamics, and each 
method has advantages and disadvantages. 

Singularity Methods 

Singularity methods typically use linear superposi- 
tions of solutions of Laplace’s equation (i.e., source, 
sink, doublet, vortex elements) to model systems that 
may include the fuselage, the fuselage wake, the rotor 
blades, and the rotor wake. Johnson 2 provides an ex- 
tensive discussion of singularity methods used for ro- 
torcraft analyses up through the year 1986. Boyd 3 dis- 
cusses other examples of analyses along these lines that 
have been published since that time. These analyses have 
shown varying degrees of success. It is apparent from 
these references that (1) one of the primary advantages 
of these methods is that they are typically computation- 
ally efficient and (2) one of the primary disadvantages is 
the inability to adequately account for viscous effects. 

CFD Methods 

In recent years, CFD methods, including methods to 
solve the full potential equation, the Euler equations, and 
the Navier-Stokes equations, have become available. 4 ” 14 
In general, the full potential and Euler methods, like the 
singularity methods, have the advantage that they are rel- 
atively efficient computationally and are quite useful in 
some applications where viscous effects are not domi- 


nant. However, a disadvantage is that, for computing 
rotor-fuselage interactional effects that include viscous 
effects, a boundary layer coupling model must be em- 
ployed with these methods. To fully integrate the vis- 
cous computation, Navier-Stokes methods should be em- 
ployed. Only a few examples of Navier-Stokes com- 
putations are present in the literature. In one of these, 
Meakin 14 used the Navier-Stokes equations to compute 
the time-accurate flowfield around a V-22 tiltrotor vehi- 
cle, including the rotor. This computation was primar- 
ily geared toward demonstrating moving, chimera grid 
technology and is not currently a practical capability due 
to the large CPU times required. In general, solutions 
to the Navier-Stokes equations for interactional aerody- 
namics problems, where everything is modeled in one 
large computation, are not currently practical for routine 
use. 

Hybrid Methods 

With the expense of Navier-Stokes methods for com- 
pete rotorcraft out of reach for routine computations, 
a practical, engineering solution is to use a hybrid ap- 
proach. In hybrid approaches, several different methods 
complement each other. For example, Steinhoff, et al 15 
combined a vorticity capturing method with a Navier- 
Stokes method to reduce artificial dissipation effects on 
rotor wake vortices, which in turn relaxes the grid reso- 
lution requirement to resolve and maintain a rotor wake 
vortex in the solution procedure. Boyd and Barnwell 16 
first introduced a hybrid method that loosely couples a 
Generalized Dynamic Wake Theory 17 ” 20 (GDWT) with 
a Navier-Stokes method. Boyd 3 extended that method 
to include both a fuselage and a rotor and computed un- 
steady fuselage surface pressures and unsteady inflow for 
a complete configuration. 

The current work uses the method of Boyd 3 and 
presents results using that method. Below, a brief de- 
scription of the method is provided for completeness. 

Computational Method 

The current computation method is a hybrid 
method that loosely couples the GDWT to a Navier- 
Stokes method, OVERFLOW . The details of this 
coupling can be found in Boyd, 3 but a brief outline is 
presented here. 

As discussed earlier, determination of the gross load- 
ing and rotor trim requires many revolutions of the ro- 
tor. As such, this computationally expensive portion of 
the method is separated from the CFD portion of the 
computation. This separation greatly reduces the time 
spent on time-accurate computations in the CFD portion 
of method. Based on the above assumption, the cur- 
rent method splits the interactional aerodynamics prob- 
lem into three distinct pieces: (1) the Rotor Loading 
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Figure 2. Current hybrid method. 


Model, (2) the Rotor/Fuselage Flowfield Model, and (3) 
the Coupling Model. The arrangement of these pieces is 
shown in figure 2. 

Rotor Loadine Model 

To reduce the computational expense of the entire pro- 
cess, a model is used to determine the loading distribu- 
tion on and the trim state of the helicopter rotor. The 
model used here is based on the GDWT as discussed 
above. This model uses a solution of the Laplace equa- 
tion for a isolated, circular wing developed by Kinner. 21 
Essentially, Kinner’s solution provides admissible accel- 
eration potential functions on the circular wing. To deter- 
mine the unknown coefficients in Kinner's solution, Pe- 
ters, Boyd, He, 18 He, 17 and Peters and He, 19 used the lin- 
earized Euler equations, the continuity equation, and spe- 
cial rotor boundary conditions, to relate the Kinner accel- 
eration potential to the induced inflow at the rotor disk. 
For the current research, the resulting closed form ma- 
trix equations are iteratively solved in conjunction with a 
modified Newton-Raphson trim technique to determine 
the unsteady induced inflow, the trim state, and the un- 
steady loading distribution of the isolated rotor. 

With the solution of the GDWT for the isolated rotor, 
the “Rotor Loading Model” portion of figure 2 is com- 
plete. In figure 2 it can be seen that the pressure (loading) 
distribution from the Rotor Loading Model is used in the 
“Rotor/Fuselage Flowfield Model”. 

Rotor/Fuselage Flowfield Model 

Now, with a known pressure distribution on the rotor 
disk, a Rotor/Fuselage Flowfield Model is used to solve 
the Navier-Stokes equations. For this model, a thin-layer, 
Navier-Stokes code (OVERFLOW 22 ) has been modified 
to include an unsteady boundary condition. For this new 
boundary condition, a cylindrical, non-rotating grid is 


used to represent the rotor. The predetermined pressure 
distribution is applied as an additional term in the energy 
equation as follows: 


A(pe 0 ) = 


A(r)AP 

Y-l 


( 1 ) 


where equation (I) is in terms of the non-dimensional 
quantities used in OVERFLOW and A(r) is the ratio be- 
tween the local actual blade area and the local compu- 
tational cell area at a given radial station on the blade. 
This ratio is used to maintain the correct overall thrust. 
The additional conservative energy term in equation (1) 
is then split into two parts. One half of the term is applied 
to the “upper rotor plane” (see figure 3b) and the negative 
of the other half of the term is applied to the “lower ro- 
tor plane”. This procedure effectively creates a pressure 
jump between two planes in the rotor grid, separated by 
an “iblanked plane” which ensures that the artificial dis- 
sipation terms, which operate on a pressure discontinu- 
ity, do not modify the input pressure distribution at the 
rotor plane. All remaining flow quantities on the upper 
and lower rotor planes are determined by averaging the 
quantities at planes “A” and “B” in figure 3b. Figure 3a 
shows a top view of the rotor grid used in figure 3b. In 
this top view, a rectangular section is used to represent 
the actual blade area, and a shaded wedge represents the 
computational area (these areas are not to scale). Only 
one blade is represented in this figure. 

For a multibladed rotor, one of these computational 
wedges exists for each blade. A radially varying, addi- 
tional conservative energy term is applied along each of 
these computational wedges for each blade. At each time 
step in the time-accurate solution procedure, the pres- 
sure jump “travels” around the rotor azimuth direction, 
one grid line per time step. This unsteady boundary con- 
dition effectively represents the rotor blades as a pres- 
sure jump traveling around the rotor azimuth on a non- 
rotating, cylindrical grid. 

Using the chimera grid techniques available in OVER- 
FLOW, the above rotor grid is combined with other grids 
which represent the fuselage and the remaining flowfield. 
OVERFLOW then solves the time-accurate, thin-layer 
Navier-Stokes equations on this set of grids, along with 
the unsteady, pressure jump boundary condition. The so- 
lution procedure is executed until the initial transients are 
removed and a periodic flowfield is obtained. 

Since the specified pressure jump was originally de- 
termined by an isolated rotor model, the pressure jump 
boundary condition does not represent the combined 
rotor-fuselage system. Therefore, once a periodic solu- 
tion has been obtained with the original pressure jump 
boundary condition, an “Inflow Correction” method is 
used to account for the presence of the fuselage in the 
Rotor Loading Model. Discussion of this method is be- 
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Figure 3. Schematic of new boundary condition. 

yond the scope of this paper, but is discussed in detail 
in Boyd. 3 Figure 2 shows the location of the “Coupling 
Model (Inflow Corrections)” portion of the model. 

With these inflow corrections, the GDWT model is re- 
executed to obtain a new unsteady pressure jump bound- 
ary condition that has been corrected to account for the 
presence of the fuselage. This cycle is repeated until 
there is no significant solution change between iterations. 

Results 

Experimental Setup 

Results from the computational method discussed 
above will be compared to experimental data. The ex- 
periments used here are discussed in other references, 3,23 
but are discussed briefly here for completeness. There 
are two experiments that are used here. The first ex- 
periment (“Experiment 1”), reported by Elliott, Althoff, 
and Sailey, 23 used a Laser Velocimetry (LV) system to 
measure the induced inflow in a plane that was one ro- 
tor blade chord above the rotor tip path plane. These 
measurements were carried out for the combination of a 
generic helicopter fuselage (known as the ROtor Body 
Interaction (ROBIN) fuselage) and a four-bladed, rect- 
angular rotor system in the NASA Langley Research 
Center 14- by 22-Foot Subsonic Tunnel (see figure 4). 



Figure 4. Laser velocimeter experiment, NASA Langley 
Research Center 14- by 22-Foot Subsonic Tunnel. 

The LV measurements were processed at an azimuthal 
resolution of approximately 2.8° . Comparisons to both 
the time-averaged and the time dependent measured data 
will be made subsequently. 



Figure 5. Unsteady surface pressure experiment, NASA 
Langley Research Center 14- by 22-Foot Subsonic Tun- 
nel. 

The second experiment (“Experiment 2”) used here 
was carried out by the third author and her colleagues, 
again using the ROBIN fuselage with the same rectan- 
gular rotor system. The primary difference in the con- 
figuration between the first and second experiments is 
that, in the first experiment, the rotor drive system was 
contained inside the fuselage shell, whereas, in the sec- 
ond experiment, the rotor and fuselage were mounted on 
separate systems. That is, in the second experiment, the 
rotor drive system was mounted to the tunnel ceiling and 
the fuselage was sting mounted on a post attached to the 
tunnel floor (see figure 5). This experiment was con- 
ducted in two phases: (1) an isolated rotor configuration 
(with the fuselage lowered to the tunnel floor) and (2) a 
rotor/fuselage configuration (with the fuselage in place). 
In the first phase of this test, unsteady inflow measure- 
ments were made at a limited number of locations on the 
advancing side of the rotor, one chord above the tip path 
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plane. In the second phase of the experiment, unsteady 
surface pressures were measured along the top centerline 
of the fuselage and at several locations on the sides of the 
fuselage. The data presented here are a small subset of 
the total data taken in the second experiment. Subsequent 
comparisons will be made to these unsteady inflow and 
unsteady surface pressure data. Table 1 lists several of 
the operating conditions and rotor parameters associated 
with both Experiments 1 and 2. 


Table 1. Operating conditions and rotor parameters. 


Property 

Value 

Blade planform 

Rectangular 

radius 

0.8606 meters 

root chord 

0.0660 meters 

tip chord 

0.0660 meters 

number of blades 

4 

root cutout location 

0.24R 

flap/lag hinge location 

0.06R 

airfoil section 

NACA 0012 

twist 

-8" 

nominal thrust coefficient 

0.0065 

solidity 

0.0977 

nominal hover M tip 

0.55 

approx, mean coning angle 

1° 

shaft tilt 

3° nose down 


Induced Inflow Comparisons 

Time- Averaged Induced Inflow 

Once the iteration procedure has concluded, as dis- 
cussed in Boyd, 3 comparisons between a number of 
quantities are possible. For these comparisons, the cur- 
rent model was executed with and without a fuselage in 
the solution procedure. As shown in Boyd and Barn- 
well 16 and Boyd, 3 the current model is also applicable to 
an isolated rotor configuration ( Le ., no fuselage). 

First, a comparison is presented between the measured 
and predicted, time-averaged induced inflow. Inflow ra- 
tio is defined as the local velocity divided by the rotor tip 
speed. The measurement data are from Experiment 1 at 
a plane that is one blade chord above the tip path plane 
of the rotor at a rotor advance ratio of jj = 0.23. The 
predicted results are from the same location above the 
rotor tip path plane and are at the same operating condi- 
tion used in Experiment 1. The rotor tip speed is used to 
make the data and predicted results nondimensional. 

Figure 6a shows the measured, time-averaged induced 
inflow ratio parallel to the rotor tip path plane from Ex- 
periment 1. These experimental data show an induced in- 
flow pattern that is not symmetric between the advancing 
and retreating sides of the rotor. For example, the contour 


line with a value of 0.015 shows that the induced inflow 
is asymmetric about the fore-aft plane of the rotor. Figure 
6b shows the predicted, time-averaged induced inflow ra- 
tio parallel to the rotor tip path plane for the isolated ro- 
tor configuration. Although the magnitudes are similar 
to the measured values, the inflow distribution does not 
match the measured distribution well. Here, unlike the 
measured data, the predicted induced inflow is somewhat 
symmetric between the advancing and retreating sides of 
the rotor. Figure 6c shows the predicted, time-averaged 
induced inflow ratio parallel to the rotor tip path plane for 
the full rotor-fuselage configuration. It is seen that the 
fuselage has a large impact on the inflow distribution. As 
with the isolated rotor configuration, the magnitude of 
the inflow matches the measured data well. In addition, 
the distribution of inflow now matches the experimental 
data well, including the asymmetric pattern seen in the 
measured data. Figure 6d shows the difference between 
the full rotor-fuselage configuration and the isolated ro- 
tor configuration. This difference plot shows the effect 
of the fuselage on the in-plane induced inflow. As would 
be expected for a fuselage, there is a deceleration of the 
flow over the forward portion of the rotor disk due to the 
upward slope of the nose of the fuselage and a subse- 
quent re-direction of the flow. Over the rear portion of 
the rotor disk, there is an acceleration of the flow due to 
the downward slope of the rear portion of the pylon. 

Figure 7a shows the measured, time-averaged induced 
inflow ratio perpendicular to the rotor tip path plane from 
Experiment 1. These measured data show several typi- 
cal features of time-averaged induced inflow. First, there 
is an upwash on the forward portion of the rotor disk. 
Second, there is an increased downward inflow toward 
the rear portion of the disk with concentrations in the 
first and fourth rotor quadrants. Figure 7b shows the pre- 
dicted, time-averaged induced inflow ratio perpendicular 
to the rotor tip path plane for the isolated rotor config- 
uration. This configuration exhibits many of the same 
features as the measured data. For example, there is an 
upwash on the forward portion of the rotor disk, but that 
upwash is not as prominent as in the measured data. Fig- 
ure 7c shows the predicted, time-averaged induced in- 
flow ratio perpendicular to the rotor tip path plane for the 
full rotor-fuselage configuration. The magnitude as well 
as the inflow distribution is well matched when the fuse- 
lage is included in the computation. Figure 7d shows the 
difference between the full rotor-fuselage configuration 
and the isolated rotor configuration. Again, this figure 
displays features that are expected due to the presence of 
a fuselage. For example, there is an increased upwash 
over the forward portion of the disk as the flow is de- 
flected upward over the nose of the fuselage, and there is 
an increased downwash at the rear of the rotor disk, just 
aft of the pylon, as the flow accelerates downward just 
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behind the fuselage pylon. 

Time- Accurate Induced Inflow 

The previous section showed that the time-averaged 
induced inflow in the parallel and perpendicular direc- 
tions (relative to the rotor tip path plane) are well pre- 
dicted by the current unsteady method. This section will 
present comparisons of the measured and predicted un- 
steady inflow data corresponding to the same flight con- 
ditions used in Experiment 1. The measured data pre- 
sented here is from the first phase of Experiment 2 (iso- 
lated rotor configuration). 

Figure 8 shows the measured and predicted unsteady 
induced inflow ratios. These inflow ratios are at an az- 
imuthal location of \\f = 84 a and a blade radial location of 
r/R = 0.80. Both the isolated rotor and combined rotor- 
fuselage configuration are shown. Both components are 
well predicted, especially the inplane component. For 
this particular location, the presence of the fuselage has 
only a minor impact on the predicted unsteady induced 
inflow. Previous literature has shown 3, 16 that these in- 
duced inflow ratios are typically well predicted over the 
entire rotor disk. 

Unsteady Surface Pressure 

In Experiment 2, unsteady surface pressure measure- 
ments were made for the same flight configuration and 
the same flight conditions as in Experiment 1. These 
measurements were made along the top centerline of the 
fuselage and at several locations on the sides of the fuse- 
lage. Comparisons are made here between the measured 
and predicted unsteady surface pressures along the top 
centerline and at several locations on the advancing and 
retreating sides of the fuselage. These pressure taps on 
the sides of the fuselage were located at several vertical 
locations and at a constant 44% of the fuselage length. 

For these comparisons, a modified pressure coefficient 
is used. This modified pressure coefficient is defined in 
equation (2) and is used to avoid numerical problems as- 
sociated with the definition of the standard pressure coef- 
ficient when the freestream velocity approaches zero (as 
would be the case in hover). 


c ' = 100 (/>-/>-) 
p ip ^) 2 


( 2 ) 


In equation (2), P is the local pressure, P U is the 
freestream pressure, p is the freestream density, OR is the 
rotor tip speed, and the factor of 100 is included for nu- 
merical convenience. For reference, equation (3) shows 
the relation between the standard pressure coefficient and 
the modified pressure coefficient used here. 


c' p = lOOpiCp (3) 

In equation (3), ^ is the standard rotor advance ratio and 
C p is defined in the usual way. 


Figure 9 shows a comparison of the unsteady compo- 
nent of the measured and predicted modified pressure co- 
efficient on the top centerline of the fuselage at various 
stations along the length of the 2 meter long fuselage. 
The location of the reference blade is plotted along the 
horizontal axis, and the negative of the modified pres- 
sure coefficient is plotted along the vertical axis. Since 
this is a four-bladed rotor, a dominant pressure pulse can 
be seen at a frequency of four pulses per rotor revolution. 
This is indicative of the four blades individually passing 
over each measurement location. It can be seen that the 
phase of each of the predictions matches the measured 
phase well; however, the amplitudes are slightly over- 
predicted. 

Figure 10 shows a comparison of the unsteady com- 
ponent of the measured and predicted modified pressure 
coefficient on the left and right sides (retreating and ad- 
vancing sides, respectively) of the fuselage at a constant 
downstream location of x — 0.8809 meters (x/L ^ 0.44) 
for several vertical locations. Again, the reference blade 
location is on the horizontal axis, and the negative of the 
modified pressure coefficient is on the vertical axis. The 
retreating side comparisons show that the unsteady pres- 
sures are slightly overpredicted, while the advancing side 
unsteady pressures are well matched in magnitude and 
phase. 


Conclusions 

A novel computational model for unsteady rotorcraft 
interactional aerodynamics has been presented. This 
new hybrid model couples a rotor loading model and a 
rotor/fuselage flowfield model in a manner that is effi- 
cient and capable of predicting time-averaged and time- 
accurate rotor inflow ratios and unsteady surface pres- 
sures on the fuselage due to blade passages. 
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Figure 6. Measured and predicted time averaged parallel induced inflow ratio from time accurate computations. 


(a) Data (b) Prediction for Isolated Rotor 







Figure 7. Measured and predicted time averaged perpendicular induced inflow ratio from time accurate computations. 
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Figure 8. Measured and predicted induced inflow in two directions for an isolated rotor and a rotor/fuselage combina- 
tion. 
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Figure 9. Measured and predicted unsteady modified pressure coefficient on the top centerline of the ROBIN fuselage, 
"x” denotes the distance in meters from the nose of the 2 meter long fuselage. 
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Figure 10. Measured and predicted unsteady modified pressure coefficient on the retreating and advancing sides of the 
ROBIN fuselage, “z” denotes the distance measured in meters from the horizontal reference line of the fuselage. 
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